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We develop a robust and scalable strongly-coupled multiscale atomistic method with electric field- 
matter interactions that can scale to multimillions of atoms. Mathematically, we need to obtain 
efficient and accurate numerical methods to compute forces and energies in large discrete systems 
where the particle interactions involve a combination of nonlinearity and long-range. This is the 
mathematical idealization of a crystalline system that is close to periodic in much of the domain, but is 
defected in certain regions. The defects have an important role in functional materials, and we are 
working to understand the structure and response of defects starting from physics-based models. 


The key challenges are 


* very large number of interactions — multi-million or much larger system size in terms of 
degrees of freedom, 


* both short-range and long-range terms, 
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* long-range converges very slowly (“every charge interacts with every other charge”), but cannot 
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which are not coarse-grained 


On the theoretical side, we have developed new continuum limits for systems with long-range 
interactions; this extends work by R. D. James (Univ. Minnesota) and co-workers. This enables the 
coarse-graining of these interactions away from defects. Future work will extend these to include 
positional disorder (due to thermal fluctuations). 
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Furthermore, while density functional theory (DFT) can’t handle charges — rigorous results show that 
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predictions. This was not part of our focus, but is the focus of much research activity worldwide for a 
variety of materials. Our methods and code will benefit from these advances. 


Finally, extremely large shared-memory systems (e.g. Blacklight and Greenfield at Pittsburgh 
Supercomputing Center, www.psc.edu) enabled our algorithms that exploit this hardware. 


2. Summary of the most important results 


One code that we have developed is a high-performance multiscale atomistics code for structural 
materials with multiple species. This builds on previous code developed by Dr. J. Knap (ARL) for 
single species systems. The capabilities of this code are unprecedented; it is now being applied by 
former PhD student Dr. Jason Marshall (now at Caltech) together with Dr. Knap to examine 
dislocation-void interaction in NiAl, a problem of interest to the Army. This is the largest, and only / 
first, calculation of this process that takes into account the full geometrical complexity and a 
sufficiently large domain that spurious boundary/image effects can be safely ruled out. It is providing 
important new physical insights that can impact materials-processing strategies. Two papers are in 
preparation on this work: one paper describing the methodological advances, and a second showing the 
application. Drs. Marshall and Knap are the authors of these papers. 
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surface with complex and spatially-heterogeneous electric fields. It has revealed unexpected new 
physics in the distortion around surface electrodes. 


The method makes certain key assumptions, namely, that the lattice is crystalline far from the defect. It 
is distorted, but affinely so, thereby preserving crystalline order, with the distortion varying on a 
lengthscale much larger than the atomic spacing. This is the key fact that enables coarse-graining in 
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bring statistical mechanics into play. 
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Abstract 


Modern technology is continually evolving and pushing the limits of our under- 
standing of materials. This is especially true at the smallest length scales, where 
phenomenological and experimental data is challenging to acquire. Computational 
methodologies are becoming more important in understanding the physics of ma- 
terials at these small length scales. Multiscale methods, in particular, are playing 
a prominent role in understanding the importance of defects and how they effect 
properties. We focus specifically on the Quasicontinuum Method (QC). We extend 
the local Cauchy-Born QC to long-range Coulombic interactions. We then use the 
method to simulate electrical and mechanical loading in an ionic solid. Next we ex- 
tend an existing non-local QC force method that was initially intended for crystals 
with a single atomic species to a multi-lattice system. The problem is treated as 
the union of many simple lattice problems and implemented in an existing compu- 
tational framework. The method is used to investigate void and vacancy defects in 
NiAl. Lastly, we develop a method for including electrostatic interactions in non- 
local QC based on the Fast Multipole Method (FMM). We use inspiration from the 
FMM to derive error bounds on our computed electrical quantities. The numerical 
error is investigated against a variety of simulation parameters and discussed. The 
final implementation is written for high performance computing applications and can 
be used to simulate a variety of materials. Overall, the newly developed methodology 
and code allows the simulation of defects at the atomistic length scale in ionic solids 
with real-space boundaries. 

Keywords: multiscale modeling, electrostatics, Quasicontinuum (QC) Method, 


ionic solids, Fast Multipole Method (FMM) 
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Introduction 


Smaller and faster are two important mantras in the design of cutting edge technolo- 
gies. These mantras are especially true in regards to electronic devices. Designing 
next generation devices to meet these goals requires exploiting the inherent structure 
and properties of specific classes of materials. Electromechanical materials, which 
have a coupling between electrical and mechanical properties, are one such class. An 
important feature of these materials is near instantaneous electrical response (speed) 
at small length scales (size), in some cases even the atomistic length scale. Next gen- 
eration devices currently being designed and developed with these materials include 
random access memories, actuators, and sensors, among many others [3, 64, 45]. 
Performing experiments on electromechanical materials, at the length scales sought 
for these devices, is both difficult and expensive, making design through experi- 
ments alone challenging. Computational methods are being actively sought by many 
researchers to augment experimental data and aid in the design of these devices. 
However, numerous challenges are inherent in any computational method that seeks 
to model materials at these length scales. 


We highlight the importance of Gallium-Nitride (GaN) as an example of an impor- 
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tant electromechanical material and the challenges in simulating it. This is just one 
example, there are literally hundreds of materials equally as important. GaN is one 
of the most important and prevalently used materials in modern technologies. The 
widespread use of GaN can be attributed to its many useful properties. Applications 
that utilize GaN material properties include LED lighting[57], solar cells[70], High- 
electron-mobility transistors (HEMT) [53], and many other technologies[60, 32, 33}. 
In fact, the 2014 Nobel Prize in Physics was awarded to researchers for their work 
on creating blue LEDs with GaN[55, 26]. GaN’s material properties and array 
of applications has led to a wide breadth of theoretical research at various length 
scales[11, 36, 81]. We only list a few sample publications at different length scales, 
there are hundreds of others that could be referenced as well. Even with this in- 
tense focus on GaN, additional research is needed, especially in regards to the role 
that electrostatic interactions play in GaN. We aim to to further investigate the 
role of electrostatics interactions, especially in systems with defects, non-periodic 
geometries, and non conducting boundary conditions. 

One of the main challenges is that continuum theories break down, for all mate- 
rials at small length scales, requiring the use of atomic theories, at a minimum, in 
areas of interest. However, modeling a complete system solely with atomic theories 
quickly becomes computationally intractable in size, because of the large number of 
degrees of freedom. The mechanics community has sought to overcome this chal- 
lenge by developing multiscale methods that bridge the gap between continuum and 
atomistic length scales, some of these methods include [78, 68, 72]. The general goal 
is to provide atomistic resolution only in areas of interest and continuum resolution 
elsewhere. One of the methods that has experienced success in achieving this goal 
is the Quasicontinuum (QC) method. The QC method can be broken down into 
multiple sub-methods [68]. While there are many flavors of QC, a common feature 
is the use of interpolation functions for atomic positions to reduce the number of de- 
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grees of freedom in the system to more manageable levels. Energies or forces in the 
system are then evaluated through numerical quadrature rules and allow solutions 
to be found through minimization routines. 

One important paradigm between these methods is local versus non local. Lo- 
cal methods, like [67], utilize the Cauchy Born rule with standard finite element 
quadrature rules to evaluate the energy. This method is considered local, because 
the quadrature rules depend only on variables at that individual quadrature point. 
An issue immediately arises, however, as the Cauchy Born rule assumes a periodic 
lattice inside each element. This assumption is valid only when the element size 
is large compared to the lattice length scale. Obviously in areas of interest, (for 
example, near defects, surfaces, domain walls, external loadings, etc.) the use of 
large elements does not provide adequate solution resolution. Researchers overcame 
this difficulty by introducing non local methods [40, 19] with different quadrature 
rules, typically based at nodes. These methods achieve full atomistic resolution, 
but are computationally and numerically more difficult to implement, because the 
quadrature rules are now non local. Essentially, the quadrature points are a function 
of variables at other quadrature points (non local) in addition to their own (local). 
This distinction between local and non local is important for later sections of this 
thesis. 

While the QC method, along with other multiscale methods, is effective for the 
analysis of many types of materials, difficulties occur for electromechanical materials. 
Materials which current QC methods can model are characterized by interaction 
energies between atoms that are short ranged and can be approximated by a cutoff 
radius. Electromechanical materials, however, are long-ranged and therefore lack 
an effective cutoff radius. Because of the lack of cutoff radius, current multiscale 
methods are not valid for this important class of materials. 

Even methods designed to specifically solve electromechanical systems have issues 
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that are difficult to overcome. The more commonly used methods include Ewald 
summations [21, 23, 74], Ewald-like summations[22, 23], and Fast Multipole Methods 
(FMM)[28]. The Ewald and Ewald-like summations assume some form of periodicity 
and do not allow easy implementation of far field boundary conditions. All of these 
methods, while computationally fast in comparison with the full summation, do 
not easily allow explicit handling of all electrostatic boundary conditions. The lack 
of explicit control of boundary conditions and assumption of periodicity inhibits 
the ability to model many geometries and defects important in materials like GaN. 
These boundary conditions are pivotal in ionic solids like GaN [23, 54, 2], especially 
in the simulation of defects. The FMM is a fantastic method, however, it is not as 
prevalently used in the atomistic and multiscale modeling field as other methods. 
Its algorithms are challenging and particularly difficult to implement in conjunction 
with other codes, especially in a high performance computing environment (HPC). 

In Chapter 2, we formally extend the local QC method [67] to handle long range 
interactions present in electromechanical materials. A more complete introduction 
and background information for the specific problem is in Section 2.2. Our new 
method uses polarization as a multiscale mediator quantity to bridge the atomistic 
and continuum length scales. We use our method to analyze and gain insight into 
the fundamental physics of electromechanical materials. This work was done in 
collaboration with my advisor, Dr. Kaushik Dayal and is published in [49]. 

In Chapter 3, we extend the non-local cluster force QC method [40] to multi- 
ple lattices. This problem was first addressed in [42], though we develop our own 
computational implementation separate from that work. The complex QC method 
essentially treats the problem as the union of a number of simple lattice problems. 
The main difference between the complex lattice case and the simple lattice case is 
that in the complex lattice case each lattice interacts with the other. We address this 
issue by creating a new C+-+ class that computes the force interactions between any 
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two lattices in a domain. This work was done in collaboration with my mentor at the 
Army Research Laboratory, Dr. Jaroslaw Knap and is preparation to be submitted 
to IJNME [50]. 

In Chapter 4, we extend the work done in Chapter 3 to handle electrostatic 
interactions. We use a combination of the FMM [28] and Reaction Field (RF) method 
(73, 27, 6] to accomplish our goal. Essentially around a point of interest, we need to 
be able to calculate the electric field. All atoms that are near the evaluation point 
are handled explicitly. Contributions to the electric field from atoms farther away 
are approximated with the first term of the multipole expansion, the dipole term. 
We take all higher order moments as a bound on the error from the approximated 
field contributions. We derive both analytical error bounds and numerically simulate 
errors with the QC method. This work was in collaboration with both Dr. Dayal 
and Dr. Knap, and is in preparation to be submitted to SIAM MMS [48]. 

Additionally, in Appendices B and C we provide additional information on using 


the QC code and its computational scaling. 
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Local Atomistic-to-Continuum Multiscale Modeling 
with Long-Range Electrostatic Interactions in Ionic 


Solids 


2.1 Abstract 


We present a multiscale atomistic-to-continuum method for ionic crystals with de- 
fects. Defects often play a central role in ionic and electronic solids, not only to limit 
reliability, but more importantly to enable the functionalities that make these mate- 
rials of critical importance. Examples include solid electrolytes that conduct current 
through the motion of charged point defects, and complex oxide ferroelectrics that 
display multifunctionality through the motion of domain wall defects. Therefore, it 
is important to understand the structure of defects and their response to electrical 
and mechanical fields. A central hurdle, however, is that interactions in ionic solids 
include both short-range atomic interactions as well as long-range electrostatic inter- 
actions. Existing atomistic-to-continuum multiscale methods, such as the Quasicon- 
tinuum method, are applicable only when the atomic interactions are short-range. In 


addition, quantum mechanics simulations via density functional models are unable 


to capture key phenomena of interest in these materials. 

To address this open problem, we develop a multiscale atomistic method to coarse- 
grain the long-range electrical interactions in ionic crystals with defects. In these 
settings, the charge density is rapidly varying, but in an almost-periodic manner. The 
key idea is to use the polarization density field as a multiscale mediator that enables 
efficient coarse-graining by exploiting the almost-periodic nature of the variation. 
In regions far from the defect, where the crystal is close-to-perfect, the polarization 
field serves as a proxy which enables us to avoid accounting for the details of the 
charge variation. We combine this approach for long-range electrostatics with the 
standard Quasicontinuum method for short-range interactions to achieve an efficient 
multiscale atomistic-to-continuum method. As a side note, we examine an important 
issue which is critical to our method: namely, the dependence of the computed 
polarization field on the choice of unit cell. Potentially, this is fatal to our coarse- 
graining scheme; however, we show that consistently accounting for boundary charges 
leaves the continuum electrostatic fields invariant to choice of unit cell. 

Keywords: electromechanics, multiscale modeling, atomistics, long-range inter- 


actions, Quasicontinuum method 
2.2 Introduction 


Ionic crystals such as solid electrolytes and complex oxides are central to modern 
technologies for energy storage, sensing, actuation, and other functional applica- 
tions. Atomic-scale defects often play a central role in these materials, not only to 
limit reliability, but more importantly to enable the functionalities that make these 
materials of critical importance. E.g., in solid electrolytes, conduction is mediated 
by charged point defects [37]; and in complex oxide ferroelectrics, functionality is 
mediated by planar domain wall defects, and loss of functionality often occurs when 
domain walls are “pinned” by charged oxygen vacancy point defects [64]. A fun- 
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damental understanding of these materials therefore requires an accounting of the 
atomic-level structure of the defects. This poses a multiscale problem: atomic-level 
resolution is required at the defect, while complex geometries and boundary condi- 
tions require the modeling of a large specimen. 

While defects play a critical role in determining properties, they occupy a tiny 
volume of the lattice; an exceedingly large fraction of the crystal is close-to-perfect. 
This feature is exploited in most leading atomistic multiscale methods such as the 
Quasicontinuum (QC) method [68, 51]: typically, an adaptive coarse-graining is used 
with atomic resolution in the vicinity of the defect and a coarse-grained description 
further away as the crystal tends to a perfect lattice. Another important aspect of the 
coarse-graining is the use of sampling or quadrature to efficiently evaluate the energy 
in the coarse-grained region. It is essential that the atomic interactions are short- 
range to allow the evaluation of the energy at the quadrature points to be efficient. 
Therefore, existing multiscale methods cannot handle long-range electrostatic /ionic 
interactions which decay as 1/r, where r is the separation between charges. 

A symptom of this difficulty with electrostatic interactions can be observed in 
standard proofs of the Cauchy-Born (CB) theorem which require the interactions to 
decay faster than 1/r? [24, 10]. Roughly, this implies that the standard CB theorem 
requires that the charge distribution in each unit cell of the lattice should not have 
net charge or net dipole character, but only higher-order multipoles. As shown in 
[35], it is not possible to define a meaningful energy density W/(-) in such a setting, 


i.e., the standard decomposition of the energy used in elasticity 
= | W(-) dQ — i boundary working 
Q aQ 


is not valid. When long-range electrostatic forces are involved, W(-) does not depend 
solely on the local value of a field (the strain field in elasticity or the polarization 
field in electrostatics). Rather, it depends on the electrostatic fields in a nonlocal 
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manner as well as boundary conditions. While not always stated explicitly, some 
notion of the CB theorem is inherent in most atomistic multiscale formulations. 

The restrictions described above on the nature of atomic-level interactions for 
conventional multiscale methods exclude an extremely large class of materials; essen- 
tially, all dielectrics, polarizable solids, and ionic solids. In dielectrics and polarizable 
solids, the non-vanishing dipole moment in the unit cell is central to the physics of 
dielectric response, spontaneous polarization. In ionic solids such as ionic conduc- 
tors, the existence of charged defects is central to enabling conduction. Therefore, it 
is essential to develop methods that can handle long-range electrostatic interactions. 

In this chapter, we present a multiscale method which is tailored to allow both 
short-range atomic interactions as well as long-range electrostatic interactions. Some 
key features of these interactions are as follows. The short-range atomic interactions 
can be highly nonlinear and involve complex multibody interactions; however, they 
are typically restricted to 2nd or 3rd nearest-neighbors in a lattice. The long-range 
electrostatic interactions have the opposite features: the interactions between charges 
are entirely pairwise, but the interactions between every pair of charges in the system 
can be non-negligible. A further important feature is that the short-range and elec- 
trostatic contributions to the total energy combine additively. These features enable 
us to leverage much of the existing work for short-range interactions. In particular, 
we can use a standard version of the Quasicontinuum method [67] for the short-range 
interactions in combination with a method that we develop for efficiently computing 
the electrostatic interactions. 

While the presentation of our method in this chapter is largely formal, key aspects 
build on — and are supported by — rigorous results of James and Miller [35] and others 
following their work, e.g. [63]. We note that seminal formal results in this topic were 
earlier obtained by [75]. While these works deal with point dipoles arranged in a 
lattice, our focus is on charges; however, due to the fact that the net charge in each 
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unit cell of the lattice is 0, many of the key results carry over largely unchanged, 
as also noticed previously by [79, 59]. The central idea which we exploit is that a 
charged lattice can be coarse-grained by introducing a polarization density field. L.e., 
electrostatic quantities which in principle require the solution of a Poisson problem 
with a rapidly and almost-periodically oscillating forcing term due to charge density 
can instead be computed with a much smoother forcing which is related to the 
polarization density field. 

Efficient and accurate methods for interactions in charged systems have a long 
history in numerical methods. The key challenge is the long-range nature of inter- 
actions: in principle, a system of N point charges requires O(N) computations to 
determine quantities of interest namely the electric field or potential. For N of the 
order required for typical problems of interest today, this is completely infeasible. 
However, the seminal and beautiful Fast Multipole Method (FMM) of Greengard 
and Rokhlin [28] provided a breakthrough in enabling this in O(N) calculations with 
a controlled error. A key strength of the FMM is that the charge distribution can 
be completely arbitrary and non-uniform; however, this generality also means that 
the method does not exploit the structure in a given problem. As noted above, in 
the problems of relevance here, the crystalline structure is lost in the vicinity of the 
defect, but large parts of the crystal are almost perfect. The FMM, however, is 
unable to exploit this structure, whereas our method is more tailored (and thereby 
also less general) and appears to scale almost independent of N asymptotically. An- 
other leading method for atomic-level calculations with electrostatic interactions is 
the Ewald method (described in e.g. [76]). It is restricted to periodic settings and 
therefore inapplicable to multiscale calculations with defects and complex geometries 
and boundary conditions. 

As mentioned above, our coarse-graining of electrostatic interactions is based on 
the notion of a polarization density field. However, it is well-known, e.g. [61], that 
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the polarization density of a periodic solid depends on the choice of unit cell. At first 
sight, this is a disturbing observation and much work has been done in the materials 
physics community on using quantum mechanical notions such as the Berry phase 
to obtain a unique choice of polarization for a given periodic charge distribution 
[61]. However, an important aspect of that approach is the insistence on starting 
from an infinite periodic solid. In both the formal calculations presented here, and 
the related rigorous calculations in the references above, the starting point is a finite 
periodic solid whose limit behavior is studied. From this “real-space” perspective, the 
boundaries of the crystal lattice enter naturally into the problem, in sharp contrast to 
starting from the infinite periodic solid where boundaries are ill-defined. The critical 
importance of the boundaries is that they, roughly speaking, compensate for the 
choice of unit cell. Ie., while different unit cell choices lead to different expressions 
for the polarization density, these also lead to different bound surface charges on the 
boundaries. When both the bulk bound charge and the surface bound charge are 
consistently accounted for in the calculations, the electric field and other quantities 
of relevance to the energy do not depend on the choice of unit cell up to an error 
that scales with size of the lattice and vanishes in the limit. Therefore, we take the 
view that a unique choice of unit cell to compute polarization density is unnecessary 
and any choice of unit cell is — in principle — equally valid. I.e., the polarization is an 
intermediate coarse-grained quantity, but there is no fundamental physical reason to 
have a specific choice. In practice, notions of crystal symmetry typically are most 
useful in selecting a unit cell. Heuristically, this perspective is comparable to the 
universally-accepted view in continuum mechanics that any reference configuration 
is — in principle — equally valid. While certain choices of reference configuration can 
lead to conceptual and algebraic simplifications, there is no fundamental physical 
preference for any specific choice. What is more relevant is that it is possible to go 
between different choices with appropriate transformations to the kinematic variables 
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and the energy densities. 


The chapter is organized as follows. 


e In Section 2.3, we formulate the problem at the atomic level and briefly describe 


the well-developed QC approach to handle short-range interactions. 


e In Section 2.4, we describe our treatment of the long-range interactions. For- 
mally, we show the appearance of the polarization as an intermediary multiscale 
quantity to link atomistic charge distributions with coarse-grained fields. We 
also examine the issue of the choice of unit cell for polarization; accounting 
consistently for the boundaries for a given unit cell does not affect the coarse- 


graining. 


e In Section 2.5, we outline the kinematic coarse-graining which follows the com- 


plex local QC method [67] and other aspects of the numerical implementation. 


e In Section 2.6, we outline the model material and its response to a variety of 


electrical and mechanical loadings. 


e In Section 2.7, we discuss various aspects of the work including open problems 


for ongoing and future work. 
2.2.1 Notation 


Throughout the chapter, bold lowercase and uppercase letters denote vectors and 
tensors. The summation convention is not used in this chapter. Sums will be explic- 
itly written out to avoid confusion except where stated. 

We define L as a Bravais lattice with three independent lattice vectors that make 


up a unit cell. 


L(e;,0) = (- ER’, 2= Sve; where v' € Z, i= 123) (2.1) 


v 


IZ 


Ny and 


Total energy 

Electrostatic energy 

Interatomic potential energy (i.e. Lennard-Jones) 
Short-range strain energy density (related to U) 
Continuum body in the current configuration 

Continuum body in the reference configuration 

Charge of the atomic species indexed by s 

Charge density field in 2 

dielectric constant for vacuum 

Dipole-dipole interaction electrostatic kernel 

Displacement field 

Slow variable representing position in current configuration 
Fast variable representing position in current configuration 
Position in reference configuration 

Intra-unit cell position of species s, defined in the reference 
Polarization density field 

Electric potential 

Electric field 

Deformation gradient 

det F 

ith unit cell 

ith partial unit cell 

Continuum material point lengthscale 

Atomic lengthscale 

Lengthscale over which continuum fields vary 

Ball of radius € 

2D disk of radius € 

Gradient with respect to the reference configuration 
Gradient with respect to the current configuration 

Surface charge due to non-neutral partial unit cells. 


represent the decomposition of (2 into the partial unit cells on the 


boundary (4), and the remainder 29 := Q\Qx; see Figure 2.4. 


2.3. Problem Formulation 


We consider a crystal occupying a region §2, composed of charged species indexed by 
s, each carrying a fixed charge Q*. The notation of species is used broadly; it refers 
to ions and electrons, as well as electron “shells” as used in core-shell models [56]. 


We assume that the charges are all point charges, e.g. nuclei, or that they can be 
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represented through a center of charge as in electron shells. Therefore, the charge 
distribution p(a) is a collection of Dirac charges. 


The total energy in the body can be written in the form below: 


Etotal = 2 Ui {riz}) 4 ; S- — (22) 


: id At €o|Ti;| 
iAj wAj 
e——_)F—- 
short-range long-range 


rij is the vector between charges i and j, and @; is the charge carried by ith atom 
of species s. The function U is the given short-range interatomic potential and can 
typically involve multibody interactions. 

We restrict our attention to zero temperature. Our goal is to find local minimiz- 
ers of Ejorqy to obtain the equilibrium structure subject to applied mechanical and 
electrostatic loadings. Brute force minimization is infeasible even for the short-range 
contributions for realistically large systems. This motivated the QC method and 
related approaches [68, 51] for short-range interactions. We will use the so-called 
local QC multi-lattice method for the short-range energy largely following [67]. As 
noted above, there are two ingredients to this multiscale approach: first, a kinematic 
condensation of the degrees of freedom using interpolations, and second, efficient 
calculation of the energy sum by using sampling or quadrature in relatively uniform 
regions. The term local refers to the fact that we use a sampling approximation 
everywhere in the specimen including at the defect core where it is likely to be quite 
inaccurate. 


We begin with the species in the reference configuration arranged in a periodic 


multi-lattice. The unit cell is denoted UO and the atomic length scale / (Figure 


2.1). In a perfect lattice with short-range interactions, the energy converges to 
J oq W(grad,, u,¢°) dV2,- Here, W is the strain energy density, and grad, u and 


C* are the deformation gradient and the “shifts” or relative displacements between 
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lattices [10, 24]. While the expression for W is algebraically involved, it is con- 
ceptually simple and comes directly from U. In a perfect multi-lattice, there is a 
well-defined notion of energy per atom since every atom of a given species is in the 
same environment. Therefore, it is possible to define an energy density by finding 
the energy of the atoms in a unit cell and dividing the cell volume, which is precisely 
W. The energy naturally depends on the shape of the unit cell and the positions of 
the different species within it, and this information is contained in grad, u and ¢° 
respectively. The QC method replaces the sum in (Equation 2.2) by sampling the 
energy density W(grad,, u,¢°) and using appropriate weights. In the more sophis- 
ticated formulations of QC, the energy is computed without this approximation in 
highly-distorted regions such as the vicinity of the defect [40]. In the local QC, the 
approximation is used throughout the specimen, including at the defect. For further 
details, we refer the reader to recent reviews of the extensive literature on applying 


QC to materials with short-range interactions [68, 51]. 
2.4 Electrostatic Interactions 


In this section, we describe the formal calculations that enable us to efficiently ac- 
count for the long-range electrostatic interactions. A key aspect is the appearance 
of the polarization density as a multiscale mediator. Because our treatment here 
is formal, we go between charge density fields and point charges as convenient by 
assuming that our calculations are valid even when the charge density field consists 
of Dirac masses. Rigorous treatments of many of the key aspects are available in the 


literature [35] and also are the focus of our ongoing work. 
2.4.1 Why the Electrostatic Energy is Long-Range 


We construct and examine some simple examples to understand why the electrostat- 


ics is denoted “long-range”. E.g., the Lennard-Jones potential has interactions that 
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Figure 2.1: Domain 2 showing the separation of scales with sample charge distribution. 


decay as r~® and these interactions nominally extend to oo. As we see, there are 


! in electrostatics. 


important differences when the interactions decay as r— 

Consider a uniform lattice and 3 cases of charge arrangement within in each unit 
cell: (i) a charge, (ii) a pair of equal and opposite charges forming a dipole, and 
(iii) two pairs of equal and opposite charges that form a quadrupole with zero dipole 
moment (Figure 2.2). 

We first consider the lattice of charges. As a rough measure of energy density, we 
compute the energy of the charge in the chosen unit cell due to its interaction with 
all the other unit cells in the body. This is the product between the magnitude of 
the charge in the chosen unit cell with the electrostatic potential created by the rest 


of the body. The potential due to a charge at a distance r from the chosen unit cell 


scales as r~'. Further, at a distance r from the chosen unit cell, we consider a spatial 
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Figure 2.2: 3 cases of charge arrangement: (i) net charge, (ii) net dipole but no net charge, 
and (iii) net quadrupole, but no net charge and no net dipole. 


region with the shape of a spherical shell with unit thickness. This shell has volume 


2 2 


that scales as r~ and therefore roughly contains r~ charges. So the total potential 


co 61,2 
r=1 Fe 


at the chosen unit cell due to the rest of the system is 5} =r ro. 
That is, the energy density of the body is unbounded in the large-body limit. The 
physical implication of this calculation is that large clusters of unbalanced charges 
have extremely high-energy and are thus unlikely to be observed in real materials. 
Next consider the lattice of dipoles. As a rough measure of energy density, we 
compute the energy of the dipole in the chosen unit cell due to its interaction with 
all the other unit cells in the body. This is the product between the magnitude of 
the dipole in the chosen unit cell with the electrostatic field created by the rest of the 
body. The electric field due to a dipole at a distance r from the chosen unit cell scales 
as r-®. Further, the shell at a distance r again contains roughly r? dipoles. So the 
total electric field at the chosen unit cell due to the rest of the body is >, 4r? = 
eae - This sum nominally also tends to infinity. However, the issue is more subtle. 
The full expression for the electric field due to a unit dipole oriented in the direction 


6:5; —38; Ba a . . 
See n;. Certain components of the summation have 


n at a position x is E; = par 
alternating sign, and this leads to conditionally convergent sums, and in general this 


sum is at the border of convergence / divergence. The physical implication of this 


LF 


calculation is that the energy density of a large collection of dipoles is extremely 
sensitive to the precise boundary conditions that are imposed far away “at infinity”. 
Alternatively, in a finite body, this says that the energy density at a given point is 
extremely sensitive to the distribution of dipoles throughout the entire body. This 
physical implication is the reason to denote electrostatic interactions as “long-range” 
namely, it is not possible to define a meaningful energy density that depends only on 
the local state of the crystal. The energy density at a given material point instead 
requires accounting for the state of the body at every other material point as well 
as boundary conditions. As we see below, this also poses practical difficulties for 
standard multiscale algorithms such as QC. This makes these methods inapplicable 
to an extremely broad class of solids: in all dielectrics, polarizable media, and ionic 
solids, the non-vanishing dipole moment in the unit cell is central to the physics of 
dielectric response, spontaneous polarization, and other key electrical properties. 
Finally, consider the lattice of quadrupoles. As a rough measure of energy density, 
we compute the energy of the quadrupole in the chosen unit cell due to its interaction 
with all the other unit cells in the body. This is the product between the magnitude 
of the quadrupole in the chosen unit cell with the gradient of the electrostatic field 
created by the rest of the body. The electric field due to a quadrupole at a distance 
r from the chosen unit cell scales as r~°. Further, the shell at a distance r contains 
again roughly contains r? quadrupoles. So the total electric field at the chosen 
unit cell due to the rest of the body is 0°, Sr? = 3, 4. This sum converges 
rapidly. This setting corresponds to the case of metals and other systems with mobile 
electrons that allow the charge to redistribute itself to “shield” the dipole moment. 
This leads effectively to short-range interactions: though nominally the interactions 
are present at all values of r, the rapid convergence of the series allows truncation 
at finite cut-off without significant error. For this reason, among others, short-range 


potentials with interactions involving only nearest- and next-nearest-neighbors are 
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sufficiently accurate to model metallic and related systems. 
2.4.2 Existing Numerical Approaches to Compute Electrostatic Interactions 


The long-range nature of the electrical interactions described above leads to prac- 
tical hurdles in atomic multiscale computations. Leading methods to handle these 
interactions are Ewald sums and the Fast Multipole Method (FMM). The Ewald 
method [29] assumes perfect periodicity. This is appropriate only for perfect crys- 
tals. Approximating defect calculations by periodic supercells has severe artifacts 
even with purely short-range interactions, a difficulty much more pronounced when 
interactions are long-range in nature. The fast multipole method (FMM) reduces 
the problem from O(N?) to O(N), but is still extremely expensive with atomic mul- 
tiscale calculations in crystals often as large as N ~ 10?!. In addition, the ability of 
FMM to deal with arbitrary charge distributions also implies that it does not exploit 
the close-to-uniform distortion away from the defect [7]. 

For short-range interactions, multiscale atomistic methods such as the QC method 
borrow ideas of quadrature rules from FEM [68, 51] to evaluate the energy at vari- 
ous sampling/quadrature atoms and then use quadrature weights. This idea depends 
critically on the energy evaluation at the quadrature point being a fast calculation. 
This is not a fast calculation if the quadrature charges interact directly with every 
other charge in the system. Therefore, these multiscale methods are applicable only 
to materials with short-range interactions. Multiscale QC-based methods for Orbital- 
Free Density Functional Theory — an empirical simplification of Density Functional 
Theory for metallic systems — use a continuous charge density rather than discrete 
point charges, but formally the issues are the same. Roughly, a predictor solution is 
patched together from the periodic solution in each “element”, and then a corrector 
solution due to the defect is superposed. The efficiency and accuracy of this approach 


requires that the corrector solution can be coarsely resolved away from the defect 
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[25]. However, in general there is a spatially-varying dipole moment in the specimen 
and zero dipole in the free space; therefore, the periodic calculation in any element 
will have large errors because it replaces this complex environment by a charge dis- 
tribution with uniform dipole density. Consequently, the corrector can require fine 
resolution over much of the domain, except in settings such as metallic systems with 
net zero local dipole where they are currently applied. 

An alternate approach to accounting for the large number of charge-charge inter- 
actions is to rewrite the problem as the electrostatic Poisson equation. However, this 
will lead to a highly-oscillatory forcing term that fluctuates at the atomic lengthscale 
while the problem is posed over the entire specimen. Therefore, this does not solve 
the essential difficulty. While numerical homogenization approaches may be feasible 
because the forcing close-to-periodic in many regions of the sample, it is not clear 
how to obtain full resolution in the vicinity of the defect. Further, the Poisson equa- 
tion can be thought of as a nonlocal constraint that must be appended to (Equation 


2.2) in the minimization and therefore the essential non-local character remains. 
2.4.3 Coarse-Graining the Electrostatic Field Energy 


We now consider the coarse-graining of the electrostatic energy. As noted above, we 
will work with a charge density field p and assume that the coarse-graining is also 
valid for point charges by replacing p with appropriate Dirac masses. Our starting 
point is to write p following the ideas of 2-scale methods |4, 72]. I.e., we consider the 
setting where the charge density varies over 2 different lengthscales. There is a rapid 
almost-periodic variation of charge density at the lengthscale of the atomic unit cell 
(denoted J). In addition, there is a much slower variation over the characteristic 
continuum lengthscale denoted L. In the language of 2-scale methods, we can write 
the charge density field as p(x, y) with y := a/I and p periodic (with period of order 


one) in the second argument. A heuristic picture is that a specifies the location of 
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the material point, and y specifies the location within the material point (Figure 
a1). 
The electrostatic field energy can be written: 


pa f CO anary 29) 


eC |x _ x’ | 


We wish to examine the limit of the energy in the following setting. We introduce 
a lengthscale € that, roughly speaking, denotes the size of the continuum material 
point. The limit of interest is then //e > 0 and «/L > 0, orl <e< L. Essentially, 
the physical interpretation of this limit is that the atomic unit cell is much smaller 
than a material point, and a material point is much smaller than the lengthscale over 
which continuum fields vary. 

We can now rewrite E as 


E= | p(x, y)p(z', y’) 
(ly) €Be(a),(ly’)eBe(ae’) (2 + ly — x! — ly’| 


w,xv’EQD 


Cave da. (2.4) 


The notation 6,.(a) denotes a ball of radius € centered at x. We note that (ly) € 


B.(a) implies y € B.i(x). 
We break up £ into 2 parts: a local term when x2 = 2’, and a nonlocal term 
when aw # a’: 


/ 
px / p(w, u)o(e.¥') 1 dV,dVy 
2eQ? WY’ Be /i() l ly _—y | 
a 


local term: x=a’ 
x, a’, ! 
+ > / A wee vy 18 dVydVy (2.5) 
xe EC Qatar! ycB. (x) ,y’€B./1 (x) |ja + yY-—a@£ = y’| 
—————————_—E— 


nonlocal term: 24a’ 


In the limit that we will take, the nonlocal term represents the interactions between 
charges that are located at different material points x, 2x’, while the local term rep- 


resents interactions between charges at the same material point. 
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We introduce some notation for what follows. We denote by U the rescaled atomic 


unit cell with characteristic dimension and volume of order 1. The atomic unit cell 


with characteristic dimensions / is denoted by /LJ. We also use U1; and /LJ; to denote 


the i-th atomic unit cell in a lattice. 


The Local Contribution of the Electrostatic Energy 


For a fixed a, the charge p(x, y) is periodic in the second argument over U. Therefore, 


we begin by rewriting the local term in (Equation 2.5) in terms of integrals over Uj: 


3 3 / pPl& ee, y’) dV,dViy (2.6) 
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The periodicity, and the fact that «/l — oo, together imply that every term in the 
sum relating the interaction between cells i and 7’ can be mapped to an interaction 
between cells 0 and some 7”. Therefore, the local term can now be written 


3 6) [. pPle, yale, Y) WV,dVy (2.7) 


rE o.y’€Be 1 (x) ly - y’| 


The factor (¢/1)? is the number of terms in the sum that are replaced, obtained from 


dividing the volume of the ball of radius €/! by the volume of 


This has the form of a Riemann sum: with e < L, the term e€° is the volume 


measure. 


3 P(x, y)p(z, y') 
Sve ( / _ l a gt) (2.8) 
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The term in the brackets is the integrand and must be well-behaved, i.e. neither blow 
up nor go to 0, in the limit 1 < ¢. The natural scaling is that the charge density 


must scale as p(x, y) = p(X, y)/l where f is the charge density on the rescaled unit 


cell LJ with characteristic dimension 1, and 1X =a. Note that p has dimensions of 


charge per unit area. While this choice of scaling may appear arbitrary, we note that 
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it can be recognized as the classical dipole scaling from elementary electrostatics. 
That is, in constructing the notion of a point dipole, one starts with charges that are 
separated by a finite distance and then takes the limit of the charges approaching 
each other. However, this limit leads to a finite dipole moment only when the charge 
magnitude is assumed to scale inversely with separation, thereby leaving the product 
of charge and separation distance finite. It is precisely this scaling which is required 
here for a finite local electrostatic energy. In our setting, if, for example, we assumed 
a fixed charge density and allowed the lattice spacing to go to 0, charge neutrality 
would give us vanishing energy. 

Using the charge scaling described above enables us to map the calculation of the 


integrand to a unit domain and gives the final form: 


A(X, y)p(X,y’ 
we \ J yeOo,y’€B. (a) ly—y'| 


The energy in this form can be readily absorbed into standard energy densities that 


arise from applying the CB theorem to short-range interactions. This term has a 
number of different names: the Madelung energy in ionic solids [38], the Lorentz 
local field, the weak-short contribution [35]. 

As an example, we replace the charge density with a set of Dirac masses represent- 


ing point charges. The charge density in the unit cell is p(X, y) = d Jen, Q°5y. (y) 


and extended periodically. The local energy has the form: 


Feat = f (x 019% + Lp(a +n) Spl) Wn (210) 
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where p(®) := dyer, Q°Ydy,(y) is the polarization of the unit cell. 


The quantities S,.S are defined as: 


ij. 4 ] 
a) De agony 

yeh zeLi\y’ zE€L1\0 

Ba Bu 


In the specific case that we have only 2 point charges in a unit cell, S vanishes and 


the local energy can be written in terms of p exclusively. The dipole kernel K is 
defined in (Equation 2.12). 

An important point above is the presence of the limit in the definitions of S and 
S. As noted previously, the full sums used above are conditionally convergent. The 
use of a limit is equivalent to enforcing a particular order of summation; in this case, 
it corresponds to using “neutral spheres” using the terminology of Ewald summation. 
Physically, it enforces that the far-field boundary conditions are set to 0. The local 
contribution then is simply the energy of a uniform lattice of charges with vanishing 
far-field electric field; the lattice is uniform because the entire lattice is located at a 
single material point. In general, there can be a non-vanishing far-field electric field 
due to the other material points and continuum-scale boundary conditions, and this 


is introduced through the non-local contribution in the next section. 
Nonlocal Contribution of the Electrostatic Energy 


We now focus on the nonlocal term in (Equation 2.5), i.e., the interactions between 
charges at different material points. This contribution provides an energy that is 
very different from the standard local continuum energies. In particular, those ener- 
gies are developed from the CB theorem that in the limit does not have any direct 
atomic interactions between different material points. Here, we have a clear nonlocal 
character to the energy. 


We first introduce some notation regarding the multipole expansion. Consider 
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the electrostatic interaction for charges located at a +ly and a’ + ly’: 


1 1, 9a ( 1 Jiew-w) 


ja +ly — a’ —ly’| ~ |e — a O(a — 2x’) \ |a — a'| 
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The operator K is the dipole kernel. 


In the nonlocal term in (Equation 2.5), in anticipation of using the periodicity of 


p(a,y) in y when z is held fixed, we reduce the integrations to unit cells 


a [ ple WAP Y) 6 ay dv (2.13) 
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Assuming a separation of scales, i.e. € < L, the periodicity of p implies that the 
interaction between charges contained in B,.(@) and B(x’) can be replaced by 
interactions between charges in unit cells at « and a’, and then multiplying by the 


number of unit cells in B.(a) and B.j)(2’). 
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Canceling the factors of | and using the notion of Riemann sums as above with e? as 
the volume measure, we can write this as a double integral: 


| ( / plz, y)olz y’) ivy) dV_ dV nt (2.15) 
ave’ CO;eA4a! ye 
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As in the local contribution, we require the integrand in brackets above to be well- 
defined when | — 0. Recall the dipole scaling p = f/l: the integrand is therefore 


well-behaved if scales as [?. 


1 
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We substitute the multipole expansion from (Equation 2.12) and notice immedi- 
ately that the first term scales independently of / and the second term scales linearly 
in 1. These would then potentially cause the integrand to diverge as | + 0. However, 


we recall that each unit cell is charge-neutral, ie. [, f(x,y) dV, = 0. This causes 


both the first and second terms in the multipole expansion to vanish. Physically, 
this means that the energy is unbounded if every unit cell is not charge-neutral, e.g. 


recalling the example in Section 2.4.1. 


Next, we consider the term $/1?7K : (y—y’) ® (y—y’). The terms containing 
y@y and y’ ®y’ vanish from charge neutrality. The only remaining terms can be 


readily written as: 


/ K(x — 2’): i, A(X, y)y av, @ | A(X" y')y! dVy | dV_ dVz 
x,2/COQ;nA-a! ye y’e 


p(x) p(x’) 


= / p(a’) - K(x — 2')p(x) dVz, dV x" 
xx’ ECOQ;nA-a! 


where the terms containing y ® y’ and y/ ® y have been combined using symmetry. 
Consider finally the terms denoted O(I?). These will all go to 0 as 1 > 0. Phys- 


ically, these terms represent the contributions from quadrupole and higher-order 


moments of the charge distribution, ie. [, p(x, y)y ® y dVy and higher-order. We 


see that these terms vanish identically in the limit. Therefore, terms of higher-order 
than dipole do not appear in the nonlocal part of the continuum energy, recalling the 
example in Section 2.4.1. In general, all short-range forces — i.e. those that decay 
faster than dipolar interactions — do not contribute to the nonlocal term in the limit 
[24, 10]. 


Finally, a long but straightforward calculation using the divergence theorem and 
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integration-by-parts gives 


: p(x’) -K(a — a')p(a) dV_ dVnr = 
ax,a’CQ 
/ div p(a’)G(a — x’) div p(x) dV_, dV 
w,xe’ELD 
(2.17) 
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— 2 | n- p(x')G(a — x’) div p(x) dV dSz 
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where G is the standard electrostatics Greens function and K := V?G from (Equation 
2.12). Note that the condition x # a’ has not been written for brevity. An important 
conclusion from the above formula is that — div p is equivalent to a bulk charge 
density (the so-called “bound bulk charge density”) and that p-n is equivalent to 
a surface charge density (the so-called “bound surface charge density”). In that 
perspective, the formula above simply gives the energy of this composite charge 
distribution using the usual Green’s function relation between charge density and 
energy. When p is discontinuous along interior surfaces, this formula gives bound 
surface charges along these surfaces. 

While we have derived the equivalent bound charges using standard integration 
formulas after taking the limit of the Riemann sum, it is straightforward to derive 
these directly. Essentially, we manipulate the Riemann sum using the standard 
approach in the Riemann sum proof of the divergence theorem (e.g., [71]). In the 
interior, we find divp appearing and the boundaries of the infinitesimal element 
canceling with it’s neighbors. On the boundary of 092, there is no cancellation 


leading to the contribution —p-n. 
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Boundary Contributions due to Partial Unit Cells 


We consider now the role of boundaries. As we have mentioned above and will 
discuss below in detail, the value of the polarization depends on the chosen unit cell. 
Boundary contributions are essential to ensure that the coarse-grained electric fields 
and other quantities do not depend on the arbitrary choice of unit cell. There are 
two contributions: first, due to polarization terminations —p-n, and second, surface 
charges due to partial unit cells on the surface that are not charge-neutral. The 
polarization terminations have already been accounted for as shown in (Equation 
2.17). In this section, we consider the case of the surface charges due to partial 
non-charge-neutral unit cells. 

For simplicity, we do not compute the total energy which will have straightforward 
but tedious cross-terms between interior bound charges (due to div p) and surface 
charges as can be seen (Equation 2.17). Instead, we compute the electric potential 
field due to the surface charges where the calculations are more transparent. In a 
formal setting, one can readily go between these calculations. 

Consider a point « € 022. At a point a’ # a, the electric potential due to the 


charges at x is given by: 


A(z,Yy) 93 
eal [ malyog) 2.18 
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Here D,(a) is a 2D disk of radius € located at a. Therefore, D.(a) x In denotes a 
squat cylinder of height Cl oriented with axis n and cross section D.(a). The vector 
mn is the unit outward normal to (2. It is implicit in the above formula that we are 
considering charges only in the partial unit cells. 

We assume that the surface is a rational plane (see Appendix A for the definition). 
The integration above in the directions along the surface (i.e., normal to 7) can then 
be reduced to an integration over a single unit cell because of periodicity in those 
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directions!. The integration in the direction along n reduces simply to the partial 
unit cells on the boundary and is therefore independent of Cl. Following the ideas 
above for the volume contributions, we can then rewrite this as an integral over a 


unit cell by putting in the appropriate factor for the number of unit cells in the disk: 


ee 
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where we have used the notation A for partial non-neutral unit cells. 
Therefore, we require that only the term independent of / from (Equation 2.12) 
appear above. Upon taking the Riemann sum and defining the surface charge density 


o(a), this gives the expected and simple result: 


(fame) 
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While we have assumed for simplicity that « 4 x’, considering the case « = ax’ and 
examining the energy would give us a local contribution analogous to that in the 


case of the bulk. 
2.4.4 Role of Boundaries in Compensating for the Non-uniqueness of Polarization 


It is well-known, e.g. [61], that the value of p in a periodic solid depends on the 
choice of unit cell. This would appear to be a fatal difficulty in using p as a multiscale 
mediator between the atomic-scale-variation of p and continuum-scale quantities. In 
the materials physics community, quantum mechanical notions are invoked to obtain 
a unique choice of polarization for a given periodic charge distribution [61]. However, 


from the perspective of the calculations above, we are simply coarse-graining classical 


! Tt is not clear to us how to proceed without assuming that the surfaces are rational planes. 
Irrational surfaces cause severe difficulties in defining surface energies even in simpler models of 
solids [62]. 
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electrostatic interactions and there is no reason for quantum mechanics to play any 
role. As we describe in this section, the difficulties noticed by [61] are entirely due to 
their starting-point of an infinite periodic solid. This makes the notion of boundaries 
ill-defined. If instead we begin from a finite solid and take the limit of lattice spacing 
being much smaller than the size of the body (the large-body limit), we see that the 
surface charges on the boundaries play a critical role. In short, while the polarization 
is itself not uniquely-defined, the electrostatic energy that comes from accounting for 
both the polarization and the surface charge is a unique quantity. While changing the 
unit cell changes the value of the polarization density, it also changes the boundary 
charge in the partial unit cells. These compensate to give the same value for the 


electrostatic energy. 
One-Dimensional Illustrative Example 


Consider a finite body 2 = (—L,L) x (—1,1) x (—1,1) with a one-dimensional 
charge distribution p(a) = po sin(27*+) (Figure 2.3). We assume that L is an integer 
multiple of J, i.e. nl = L. Guided by the multipole expansion, we compute the dipole 
moment as the leading contributor to the behavior of the bar without using the fact 
that the charge distribution is in fact periodic in (2. We then compare this to the 
result obtained by using polarization density field that is defined on the unit cell. 


The total dipole moment of the bar P := f, p(r’)r’ dV,. evaluates to (1,0,0) x 


—ApoLl 


TT 


The polarization density p in a single unit cell [0 = (a,1+ a) x (—1,1) x (1,1) 


is defined as p := STC Sven p(r’)r’ dV, It evaluates to (1,0, 0) x =pot cos(27$). 


volume(l 


This is a classic example showing that p depends on the chosen unit cell, here 
parametrized by a [61]. From (Equation 2.17) and the associated discussion, we 


have no bulk charge because p is identical in each unit cell, but there is a surface 


2pol 


charge density given by —p-n. Therefore, there is a total charge of ““* cos(274) 
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Figure 2.3: Charge distribution in the 1D illustrative example to show the effect of bound- 
aries. 


at +L, and —2e0! cos(27%) at —L. However, there are partial unit cells at each end: 


(-L,-L +a) and (ZL —1+a,L), and these are not charge neutral. The charge 


2pol 
Tv 


in these cells evaluates to 2 (1 — cos(274)) and (1 — cos(2m4)) respectively. 
Therefore the total charge at each end, both from the partial unit cells and from 
—p-n, is 2001 These equal and opposite charges are separated by a distance 2L. 


Stele . Note that we have errors up to order 


Therefore, this is a dipole of strength 
because the charges due to the polarization terminating on the surface are separated 
by L —1; however, the key point is that in the limit of |<< L, we recover the dipole 


P. 
The General Case 


The key lesson from the example above is that it is critical to account for the charge 
in the partial unit cells on the boundary. This ensures that the coarse-graining that 
exploits the polarization as a multiscale mediator does not depend on the choice of 


unit cell. We now examine this in the 3D setting. 


First, we decompose {2 into (24, with only complete unit cells, and Qy with the 


surface layer of incomplete unit cells, Figure 2.4. 


We now consider the unit cells adjacent to (24. Our goal is to modify these 


dl 


Figure 2.4: Decompose 2 into 25 and Qy. 


unit cells in various ways, and show that the resulting changes in surface charge and 
polarization density balance each other. 

An important element of our strategy is to transform the given unit cell to a unit 
cell that has a face parallel to the surface under consideration. Appendix A shows 
that a unit cell with this property can always be found when the surface is rational. 
As in Section 2.4.3, we restrict attention to rational surfaces. 

For this special choice of unit cell, we now consider the changes in surface charge 
and polarization density for various operations. We use the notation that the lattice 
vectors tangential to the surface are hy and h3. We note that the volume of the unit 
cell can be written hy x h3- hy = |he x h3|n- hy. 

First, consider translations of the unit cell as shown in Figure 2.5. 

Consider Figure 2.5a. A translation along h, causes an increase in the un- 
compensated surface charge area density Ag = ear ‘i @P in the partial unit 


cells. The increase in the polarization density in the translated unit cell is Ap = 


PeSeaTE Ty ( te py — to py). From the periodicity of the charge distribution p(y + 


h,1) = p(y), we have Ap = hy, te p. Therefore Ap-n = Ao. 


1 
|hoxh3|n-hy 


Consider Figure 2.5b. A translation along either hz or h3 causes no change in o. 


The change in polarization is h» te p. Therefore Ap-n = 0. 


1 
|hoxha|n-hy 


In general, one can have a translation along any direction. Such a translation 
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can be decomposed into components along the lattice directions and the calculations 


above applied in succession to each direction. 


a) 


Figure 2.5: Change in charge and polarization due to a translation along h; and ho 
respectively in the special choice of unit cell. 


Second, consider distortions of the unit cell as shown in Figure 2.6. 

From reasoning following very closely the previous case of translations of the unit 
cell, we find that the relation Ap-n = Ao holds here too. 

Third, consider changes in the unit cell due to remapping of the lattice vectors 
as shown in Figure 2.7. 

As noted in Appendix A, the relation between basis vectors {f1, fo, f3} and 
{hi, he, hs} must be of the form: f; = >), uh; with p) a matrix of integers with 
determinant +1. An example of such a remapping is in Figure 2.7a. 


As shown in the example in Figure 2.6bcde, regions of the original unit cell are 


mapped to the new unit cell. For instance, maps to (1b), maps to (2b), and 


maps to (3b). Each of these regions is translated by an integer linear combina- 
tion of {hi, hz, h3}; however, each region may have a different integer combination 
translation. In addition, the uncompensated unit cell on the boundary also increases 
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Figure 2.6: Change in charge and polarization due a distortion of the unit cell in the 
special choice of unit cell. 


in extent (Figure 2.6f). As above, for a periodic charge distribution, when a charged 
region is translated by an integer multiple of a lattice vector, the consequent change 
in polarization is simply the total charge in the region times the translation distance. 


Therefore, 


1 
Ap = ? ’n, | ths | 291. 
P |hy x hgln- hy (Sum fier oe 2 re 3 Pe ( ) 


git 8 ; 
, where v;,V7,V? are integers. 


The change in the uncompensated charge is simply Ao = aa wii J4, e which 
is related to the extent of the translation along h,. This matches precisely with Ap-n. 

In the general case of modifying a given unit cell to another shape by any com- 
bination of the mechanisms studied above, we can conceptually consider mapping 
the given unit cell to the special unit cell with a surface-parallel face, conducting the 
modifications with the special unit cell, and then mapping to the desired final unit 
cell. Of course, in practice none of this need be done; as long as we are assured that 
changes in the unit cell are compensated by boundary charges appropriately, we can 
directly modify the unit cell as desired. 
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Figure 2.7: Change in charge and polarization due a remapping of the unit cell from the 
special choice of unit cell. 


In the interior of the body, div p changes by O(1) when the unit cell is changed. 
This follows directly from the definition of p in (Equation 2.16) and the chain rule, 
using the fact that X/] = a. Therefore, the bound bulk charge density is the same 
in the limit of 1/L < 1. 

We note potential connections to ideas of Null-Lagrangians in the issue of a unique 
definition of the polarization [20]. Essentially, one can have different expressions 
for the free-energies, but these lead to the same Euler-Lagrange equation but with 
different boundary conditions. Similarly, different choices for the unit cell leave 
the bulk electrostatics unchanged, because the change in the bulk polarization is 


compensated by the boundary contributions. 
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2.5 Numerical Implementation 


Our numerical implementation for the short-range interactions follows closely the 
standard QC, e.g. [67]. In standard QC, there are two essential steps: first, a re- 
duction of degrees of freedom by interpolation, typically using linear shape functions 
inspired by finite elements, with atomic resolution in critical regions and coarse- 
grained elsewhere (Figure 2.8, Figure 2.9); and second, a fast estimation of the 
energy (or derivative of the energy with respect to the retained degrees of freedom) 
using Cauchy-Born sampling. Among the many variants of QC, we use the local QC 
for multi-lattices, following [67]. Local QC refers to the use of sampling everywhere 
in the specimen, not only in the coarse-grained region but also in regions with atomic 
resolution in the interpolation. 

Dielectrics require a multi-lattice description because a dielectric response re- 
quires at least two charges in the unit cell that move independently to change the 
polarization in response to electromechanical fields. Essentially, the multi-lattice de- 
scription uses F’ to track the deformation of the unit cell, and a set of vector-valued 
fields ¢* that track the position of individual species s within the unit cell. We use 
linear interpolation for the coarse-grained displacement field uw implying that F' is 
constant in a given element. To be consistent with this spatial variation of F’, we 
use the same “constant in each element” interpolation for ¢°. If the element were 
of infinite extent, this choice of interpolation would ensure that the energy density 
converges to the standard Cauchy-Born theorem. In the local QC approximation, we 
estimate the energy of a given element by finding the energy density of atomic unit 
cells at the selected quadrature / sampling points and multiply by the appropriate 
weight. 

This interpolation of the kinematic variables F’, ¢* implies that the polarization 


p is also constant in a given element. Therefore, in terms of effective bound charges, 
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we have no bound bulk charges (divp = 0), and we have surface charge density 
(pi — pz): n at the element faces. In the coarse-grained local QC approximation, 
the evaluation of electrostatic fields consists simply of finding the fields set up by 
the charge distributions. One element of electrostatics is that the electric potential 
and field are naturally posed in the current configuration. Therefore, for numerical 
updates, we compute the electrostatic fields in the current and then pull back to 
the reference using standard electromechanical transformations [80]. We assume in 
the remainder of the chapter that we are dealing with point charges that can move 


around, but do not change their charge. 
2.5.1 Energy Minimization through Gradient Descent 


Our interest is in finding energy minimizers to the coarse-grained problem. Under the 
local QC approximation, the coarse-grained problem can be written as the standard 


continuum energy for electromechanical solids [65, 80]: 


B= ; W (grad, u, ¢°) dVng + 2 | | grad, ¢/? dVx (2.22) 
2% RS 
—————————————————— 
short-range energy long-range (non-local) energy 
div, grad, @ = div, p (2.23) 


The nonlocal expression in (Equation 2.17) can be transformed to the form above 
by first noting that the right side of (Equation 2.17) is entirely in terms of the 
electrostatic Greens function G. Therefore, the electrostatic field ¢ can be defined as 
the solution of the electrostatic equation with charges given by — div p and p-n. The 
energy density is then simply | grad, ¢|?. The surface charges due to p-n as well as 
due to incomplete unit cells appear in the boundary conditions for the electrostatic 
equation. The local contribution of the electrostatic energy has been absorbed into 
the short-range term. Note that the non-local energy integral is posed in the current 
configuration. 
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We use a gradient-flow evolution to find the (local) minimizers following [80, 84]. 
The independent variables that remain are u and ¢*. The polarization is completely 
defined in terms of unit cell geometry w and intra-cell positions of the charges ¢°, in 
turn defining the electrostatic potential ¢ through the electrostatic Poisson equation. 


Therefore, taking variations: 


grad,, u — grad,,u+negrad,,v , C° > ¢° +70" (2.24) 


Additionally, the variation in the electric field is grad, ¢ > grad, @ + 7 grad, wv but 
this variation is constrained to variations in p + p+ nq by the Poisson equation, 
i.e. div, grad, ~ = div, q. 


The definition of gradient flow results in the following statement: 


: O 
A (7 ut Ss Cs ) dn, = — ane (grad,, u + 7grad,,, v, C° + 70°) 
0 8 


n=0 


OW OW 
7 =: -O | d 
(seat a grad, v + d. aes Vie 


+ €9 i. grad, @- grad, w dV, 
R3 
(2.25) 


The nonlocal integral over R? can be transformed by multiplying div, grad, wv = 


div, q by ¢ and integrating over R®. Using integration by parts on the left side and 


then pulling it back to the reference gives: 


| grad, d- qJ dVz, = | grad, @- grad, w dVz (2.26) 
on R3 


where J is the Jacobian of the deformation. This is substituted for the nonlocal 
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integral in (Equation 2.25) to get: 


; OW 
U-v+ *+@° | dVa, = >a : nv ds, 
[, ( du : ° Jan, O(grad,, U) : 


aw 
_ i | a 
[, («iv ° A(grad,, a) eve 


+f TS 0° dVz +e f erad, d- qJ dV, 
N20 2 


(2.27) 


We now examine the relation between q and the variations v and 8. The polar- 


ization density p is defined by 


1 s s 
P= Ta 2" F¢ (2.28) 


Note that ¢° is the position of species s in the reference configuration, while p is the 
polarization density in the current configuration. Similarly, Vo is the volume of the 
unit cell in reference configuration. Hence, both Vo and ¢* are pushed forward to the 
current. 


Taking the variation of p 


p+nq = ae Q*(F + nG)(¢* + 16°) (2.29) 


1 
det (F + nG)\ 


where G := grad,, v. Noting that (det(F'+nG))~' = det(F)~' detUI+nF'G)*! = 
det(F)~' (1 — ntr(F-1G) — O(n?)), we find: 


q= ae Q* (FO° + ¢°- grad,, v — tr(F grad,, v) FC’) (2.30) 


iach 


after ignoring terms of order higher than linear in 7. Using this expression for q, we 


obtain the following expression for [ gq, grad, ¢- qJ dV, after using integration-by- 
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parts and the divergence theorem: 


1 
rad, 6: qJ dVz, = > : | rad, 6: (F'@*) dV, 
[is @-qJ dVa, % 248 Bade d (FO) Wen 


diva, (grad, 6¢*) - ¥ in) 


2 


i] 
+— rad, 0¢*° : nv dSy -f 
me (fe oe 


(2.31) 
_ s : s —T. 
V > Q (f. grad, @:-(FC*)F : nv dSx, 


-{ div,, (grad, ¢- (F'¢*) F-')-» Wa) 
2 


Defining A® := grad, 6¢°—grad, ¢-FC°F~", collecting terms in (Equation 2.31, 
Equation 2.27) and localizing using the arbitrariness of the variations gives us the 


compact form: 


:, . OW Q* 
i= divs, (as ads as oD 4’) . CS De. cov, grads @-F (2.32) 


with the boundary term: 


OW Ce _ 


These equations provide the gradient descent equations for the fields w and ¢*. The 
equation for u involves a standard mechanical stress as well as a so-called Maxwell 
or electromechanical stress. The equation for ¢* is a local equation (i.e., not a PDE), 
but is coupled through W and the electric field to the nonlocal electrostatics and the 


PDE for momentum balance. 
Electromechanical Transformations to the Reference Configuration 


For simplicity, we solve the electrostatic Poisson equation in the current configuration 
for the electric potential ¢(a) and compute the electric field grad,, ¢(a) and rewrite 
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x = x(xq) from the deformation map. In our setting, the energy (both short- and 
long-range) and electric fields are based entirely on the atomic position in the current 
configuration. Therefore, in common with much of continuum mechanics, the refer- 
ence configuration can be considered simply a change of variables for bookkeeping 
convenience. It follows that the definition of electrical quantities in the reference 
configuration are not essential. This can also be readily observed from (Equation 
2.31); the physical quantity of interest there is grad, ¢(a), and while it is certainly 
possible to define the electric field in the reference, it is not physically important. 

This has led to a number of different proposals for referential electrical quanti- 
ties in the literature. For instance, our (Equation 2.28) suggests that the reference 
polarization po(ao) is given by p(a(ao)) = det(F)~'F'po(axo), ie. the polarization 
transforms as material line elements that carry charges but with an additional factor 
accounting for volume changes. 

In [80], they assume instead p(a(ao)) = det(F)~'po(ao), and further assume 
that the reference electrostatic potential ¢9(ao) = ¢(a(ao)); this gives them that 
the electric field transforms as line elements using the identity grad, = F~7 grad,,,. 
This is consistent with the differential geometric notion of the electric field as a 
1-form, i.e. it is a quantity that is integrated along lines. 

In [58], following [18], they use yet another transformation; they work with electric 
displacement D and electric field E as primary variables rather than polarization. 
These variables are motivated by the fact that the continuum problem posed in po- 
larization and electric potential leads to a saddle-point variational problem, whereas 
the minimization structure is preserved in D and E. The electrostatic equations 
in these variables are divD = 0 and curlE = 0. For E, these imply it should 
transform as a material line element and matches with [80] as well as the differential 
geometric notion of the electric field as a 1-form. For D however, these imply a 
transformation D = det(F)~!F~? Dp following the differential geometric notion of 
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the electric displacement as a 2-form, i.e. a quantity that is integrated over surfaces; 
this is obviously identical to the standard stress transformation. This differential 
geometric notion is also implicitly exploited in Section 4 of [41] in finding the appro- 
priate averaging for the different field quantities. An essential practical advantage of 
this transformation that is exploited by [58] in their homogenization analysis is that 
the equations retain their structure in the reference configuration, i.e. div, Do = 0 
and curly, Eo = 0. 

While the transformations proposed by these other workers are physically appeal- 
ing for the reasons mentioned above, the transformation implied by the microscopic 
model of the polarization is also physically motivated. These different transforma- 
tions are not consistent with preserving the relation D = ¢)F + p between corre- 


sponding quantities in the reference. 
2.5.2 Local Quasicontinuum for Multi-lattices with Short-Range Interactions 


We use standard finite element linear shape functions, N, defined at the nodes a to 


approximate the displacements u of the lattice vectors: 


Us ». UaNa => grad, u © > Ug grad, Na (2.34) 


where u, are the nodal displacements. Defining B as the Piola-Maxwell stress tensor, 


dive, DB = diva, (as aa + co) 4") (2.35) 


we can write (Equation 2.32) as 0 = divz, B. Standard nonlinear finite element 
methods can now be used; the key difference is that we need to compute the elec- 
tromechanical contribution to the stress at every iteration. At the same time, we 
also iterate with respect to ¢* noting that our interpolation for these variables is 
piecewise constant, i.e. constant in a given element. 
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There are two short-range calculations: first, the Piola-Maxwell stress tensor, 
and second, the minimization over ¢*. Following the complex local QC method [67], 
the deformation gradient and ¢° are related to the atomic displacements through the 
Cauchy-Born rule. We assume that the energy of each element can be approximated 
by assuming a homogeneous deformation of the crystal through the deformation 


gradient. Given an interatomic potential U, we use 


OW 1 2 OU Or 


= ees 9. 
Ograd,,u  2V Or Ograd,,, U ee 


atoms in cutoff 


where r are the positions of the atoms and are obtained from uw and ¢°. Similarly, 


we can compute 


ow 1 QU Or 


acs 2V Or OCs 


atoms in cutoff 


(2.37) 


The minimization over the nodal values of wu and the element values of ¢* are con- 


ducted in a coupled manner. 
2.5.3 Coarse-graining of the Long-Range Electrostatic Interactions 


The energy minimization calculation in (Equation 2.32) requires the electric field for 
a given distribution of charge, or equivalently for a given p. As noted above, we have 
a constant F and ¢* field in each element, i.e. they are discontinuous only along 
element boundaries. Consequently, p as defined in (Equation 2.28) is also constant in 
each element and discontinuous along the element boundaries. These discontinuities 
in polarization result in surface charge densities (p; — po) -n. We compute the 
fields due to these relatively simple charge distributions using direct Greens function 


integrations. 
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2.6 Crystal Free Surface Subject to Inhomogeneous Electric Fields 


We apply the methodology described above to a simple setting of a crystal with 
a free surface subject to an inhomogeneous external electric field due to a point 
charge above the surface. We use short-range potentials based on the bi-species 
Lennard-Jones model used for Ni-Mn by [31], giving a tetragonal neutral lattice 
with a body-centered ion and a charged shell. The polarization is oriented along 
the tetragonal direction when external fields are absent, thereby providing a spon- 
taneous polarization. This provides a simple model of many widely-used perovskite 
ferroelectrics such as barium titanate and lead titanate. 

Given that this is a model material rather than numerically accurate, we aim 
to elucidate the physics of electromechanics. There are two independent ratios of 
interest. One is the strength of the ionic charges v. the strength of the external 
charge. The other is the strength of ionic interactions v. the strength of the bonded 
interactions. We explore the effect of the former by increasing external charge while 
holding ionic charges fixed. We explore the effect of the latter by scaling the elec- 
trostatic interactions. This enables us to test a class of materials that range from 
non-ionic to ionic. 

We consider a specimen with the spontaneous polarization oriented tangential 
to the surface. A single point charge is placed near the surface. In all of these 
examples, atomic resolution is used in regions of interest — particularly beneath the 
external charges — while coarser resolution is provided everywhere else. A fine mesh 
is introduced near the point charge with coarsening throughout the rest of the body. 
Figure 2.8 shows the full mesh, while Figure 2.9 shows the zoomed in atomistically 
resolved portion of the mesh with all atoms plotted in the current configuration. The 
black atoms are nodes while the green atoms are constrained by the interpolation. 


We conduct three calculations with the electrostatic interactions scaled by factors 
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of 1, 2,4 respectively, while holding the charge location and strength fixed. For the 
given charge strength, a scaling by 1 produces a surface distortion near the point 
charge, while the scaling by 2 and 4 produce nucleation-like events”. Figure 2.10 
shows the stress and polarization fields near the point charge for electrostatic scaling 
of 1. Similarly, Figure 2.11 shows the stress and polarization for the electrostatic 
scaling of 4. 


Mesh Plot of Initial QC Mesh 
300 


0% 
400 


DV NVVUPEEee a a — ee) — 
INN \ 


Y-Location 


PKERNAND 
NANAAANAAANN 
VNZAAANZAZAN 


200 400 600 800 1000 1200 1400 
X-Location 


Figure 2.8: The mesh of the entire specimen. The lengthscales are angstroms. 


The next set of calculations replaces the single point charge by multiple point 
charges (2 and 4 respectively), keeping the total charge the same as the case of the 


single point charge and with spacing between charges on the order of the height above 


2A nucleation-like event refers to localized switching of polarization. We do not expect our 
method to be qualitatively accurate beyond this regime. 
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Figure 2.9: Close up view of the atomically-resolved portion of the sample. Without 
applied fields and loads, the surface of the specimen is flat. The surface feature is caused 
by an applied electric field. 


the surface. In subsequent sets of simulations, these charges were moved closer to the 
surface until eventually a nucleation-like event occurred. Stress plots are provided for 
the 2 and 4 charges at various distances from the free surface in Figure 2.12-Figure 
2, Lo: 

An interesting feature is the presence of two stress lobes beneath the surface of 
the material. To further investigate this, we replace the point charge by a dipole 
(two charges with opposite signs). The resulting deformation as seen in Figure 2.16 
is a depression instead of the up-down pattern seen in loadings with charges of the 


same sign. Additionally, the two stress lobes previously seen disappear and merge 
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into a single stress lobe beneath the surface. 

The stress lobes in non-dipole loadings do not appear to match continuum phase- 
field calculations, Figure 2.17 following [82]. These continuum calculations are based 
on linearized electromechanics whereas the atomistic calculations in this chapter are 
inherently large-deformation. ‘The small-deformation approximation implies that 
electric fields are computed in the reference as well as potential large rotations being 
ignored. Further, the material model is also linear elastic. To examine if these 
assumptions are responsible for the inability of continuum models to capture this 
feature, we perform the following tests. First, we use precisely our approach except 
that the electric fields are computed in the reference as in linearized electromechanics. 
The resulting stress field still shows the double-lobe structure, though with a slightly 
reduced magnitude (Figure 2.18). Second, we examine the geometric linearization 
of the strains; Figure 2.19 shows a plot of the normalized difference between the 
non-linear and linearized €22 strain measures which are far too small to be the cause. 
Third, we test against a fully atomic level description without any coarse-graining. 
While the system size is small due to computational limitations of the electrostatics, 
we see a 2-lobe structure Figure 2.20. Algorithmic reasons prevent us from computing 
the stress and deformation gradient in the fully-atomic setting, so we instead examine 
the change in energy of each atom from the perfect crystal. While qualitative, this 
calculation supports the view that the double-lobe structure is an atomic effect that 
we are able to capture despite the local QC approximation. 

Regarding the lack of agreement with continuum phase-field, we conjecture that 
this may simply be due to the gradient penalty preventing the development of fine- 
scale features such as the double-lobe. Our approach can be viewed as a phase- 
field method with no gradient penalty. In typical phase-field modeling, the gradient 
penalty is taken to significantly larger than justified by domain wall widths, because 
the goal there is to predict microstructure and not defect structure. This may cause 
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fine-scale features such as the double-lobe to be washed out. Of course, it is also 
debatable whether our fully local QC model is appropriate to model such fine-scale 
features. 

Finally, we examine the effect of a mechanical indenter that is pressed into the 
surface. Plots of the stress and polarization are in Figure 2.21. These are largely as 
expected. 

We note that in these calculations, the electrostatics appears to induce a length- 
scale. However, neither classical electrostatics not the local QC — which is essentially 
classical continuum mechanics with an atomically-informed constitutive model — has 
an intrinsic lengthscale. The induced lengthscale from the electrostatics is nonlocal 
though problem-dependent, i.e. it depends on sample size, boundary conditions, and 


so on. 
2.7 Discussion 


We have presented a multiscale atomistic method for ionic solids and other ma- 
terials where electrostatic interactions are long-range. Our approach is based on 
using the polarization field as a multiscale mediator between atomic-level rapidly- 
varying charge distributions and the continuum scale electrical quantities. Our 
coarse-graining strategy relies on ideas developed by [35] and others that followed 
them [63, 80, 59]. The method that we have presented enables QC and other mul- 
tiscale methods to go beyond purely short-range interactions, that are characteristic 
of many structural materials, to functional and electronic materials of interest today. 

The method presented here is a first attempt towards the goal of understanding 
the multiscale electromechanics of defects in ionic solids. Consequently, there remain 


many important open questions. 


1. The presentation here is based on energy minimization which implies zero tem- 
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perature; recent methods based on extending this to finite temperature in the 
setting of short-range interactions can perhaps be adapted to our setting [43]. 
An alternate approach is [69], where the powerful method of Effective Hamilto- 
nians has been applied to study the finite temperature behavior of ferroelectrics 


in the local QC setting. 


. We have started from an atomic viewpoint in this chapter whereas QC methods 
have been demonstrated for the orbital-free density functional theory — a ver- 
sion of (ground state) density functional theory that is restricted to metals with 
mobile electrons. However, variants of density functional theory have difficulty 
with charged defects as rigorously demonstrated in [12]. In addition, density 
functional theory has important qualitative failings in computing bandgaps, 
van der Waals interactions, and so forth [39]. Therefore, it appears to be more 
useful to begin from the atomic level with well-calibrated and trustworthy po- 


tentials. 


. We have worked within the local QC setting which can alternately be considered 
a standard finite element approach with the constitutive relation being drawn 
from atomistics rather than a prescribed function. An important further step is 
to couple the coarse-grained model with a region with truly atomic resolution in 
the vicinity of the defect. In this regard, the key coarse-graining ideas presented 
here carry through to that setting. This open question is an area of our current 


research and we expect to report on this in the near future. 


. Our use of linear interpolation restricts us from capturing potentially impor- 
tant effects based on strain gradients, in particular the phenomenon of flexo- 
electricity that can be relevant at small scales [47]. However, this phenomenon 


is relevant at large strain gradients that typically can only occur in localized 
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regions of the sample. Therefore, a QC method that fully resolves the defect 


region can potentially capture this phenomenon even with linear interpolations. 


5. Our coarse-graining of the electrostatics in this chapter assumed complete sep- 
aration of scales. However, for a real calculation, this will naturally not be true. 
It is therefore important to study and quantify the errors when the separation 
of scales is large but not complete. This will enable the construction of an algo- 
rithm with controlled error tolerance. Directly related to this is the derivation 
of a rigorous limit as opposed to the formal presentation in this chapter. These 
related open questions are an area of our current research, following techniques 


in [63). 


The example that we have studied of the double-lobe structure beneath a point 
charge provides an important motivation for the further development of our method. 
Atomic-scale features such as the double-lobe are of importance to phenomena such 
as domain nucleation and cannot be captured by continuum models. On the other 
hand, the limited size of our admittedly-unoptimized brute force atomic calculation 
shows the need for coarse-graining efforts. 

Finally, we have examined the issue of whether polarization is unique. As noted 
above, evaluating the classical definition of polarization provides a unit cell-dependent 
quantity. This is obviously a potential disaster for a coarse-graining scheme based on 
the polarization field. The materials physics orthodoxy appeals to quantum mechan- 
ical concepts to fix a unique value of the polarization [61]. Continuum mechanicians 
have used variational notions to achieve similar ends [59]. The view that we have 
taken in this chapter is that there is simply no need to have a unique value for the 
polarization. When partial unit cells that provide boundary charges are accounted 
for consistently, the coarse-grained electric fields and other relevant quantities are 


independent of the choice of unit cell. The polarization is merely a multiscale in- 
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termediary. In continuum theories of electromechanical solids, e.g. [80, 65], the 
polarization appears in the standard local free energy density, not only in the elec- 
trostatics. In the local free energy density, the polarization can be considered as 
a quantity that tracks the atomic positions rather than relevant to electrostatics. 
In that perspective, the specific choice of unit cell is irrelevant. Broadly, the view 
advocated in this chapter is in the spirit of classical continuum mechanics: the polar- 
ization field simply provides some information about the atomic-level, and one can 
take any choice as long as consistent transformations between energy, kinematics, 
and boundaries are respected. The immediate analog in continuum mechanics is the 
freedom in the choice of reference configuration and the corresponding value of the 
deformation field and strain energy density response function as long as care is taken 


to define suitable transformations between different choices. 
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Figure 2.10: Stress and polarization due to single point charge with an electrostatic scaling 
of 1. 
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Figure 2.11: Stress and polarization due to single point charge with an electrostatic scaling 
of 4. 
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Figure 2.12: Stress due to two point charges with an electrostatic scaling of 1. 
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Figure 2.13: Stress due to two closely-spaced point charges with an electrostatic scaling 
of 1. 
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Figure 2.15: Stress due to four closely-spaced point charges with an electrostatic scaling 
of 1. 
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Figure 2.16: 
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Figure 2.17: Stress due to a point charge computed using a continuum phase-field model 


[82]. 
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Figure 2.18: Stress due to a point charge with electrostatic fields computed in the reference 


configuration. 
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Figure 2.19: Normalized Difference between Nonlinear and Linearized Strain Measure (€22 


component). 
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Figure 2.20: Fully atomic calculation for a sample with a point charge. The plot shows 
the energy difference from the state without any external charge with the electrostatic field 
energy subtracted. The entire computational domain is shown. 
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Figure 2.21: Stress and polarization due to mechanical indentation with no applied electric 
field. 
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3 


Development of a Fully Non-Local Quasicontinuum 
Method for Simulating Atomistic Defects in 
Materials with Many Atomic Species 


3.1 Abstract 


Atomistic defects (e.g. voids, vacancies, domain walls, etc.) play a vital role in 
continuum level properties in materials ranging from high strength alloys like NiAl 
and TiAl to ionic solids like Ceria and Barium Titanate. Full atomistic resolution is 
required around these defects to fully understand the underlying physics, however, 
atomistic resolution cannot be obtained everywhere because of computational lim- 
its. Multiscale methods were developed to make problems tractable by providing 
atomistic resolution near defects and coarse-grained resolution elsewhere. A variety 
of methods have been proposed and developed, but relatively few have the ability 
to simulate materials with multiple species of atoms. In this work, we extend one of 
these multiscale methods, the non-local quasicontinuum method, to multiple species 
of atoms. We extend the method by treating the problem as the union of multiple 


single species (simple) lattice problems. Treating the problem in this manner allowed 
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the leveraging of existing simulation codes in the extension of the method. We use 
this newly developed method to simulate void relaxation in NiAl. 
Keywords: non-local Quasicontinuum method; multi-lattice multiscale model- 


ing; voids; NiAl 
3.2 Introduction 


Atomistic defects play a pivotal role in the structure and functional use of contin- 
uum level material properties. For example, the motion of charge defects enables 
the functionality in next generation solid oxide fuel cells and batteries; the motion 
of domain walls in ferroelectric materials enables the storage of state used in ferro- 
electric RAM; and void defects in high-temperature alloys play a role in ductility 
and plasticity important for lightweight turbine blades. Analyzing and understand- 
ing the underlying physics of these defects is critically important in device design, 
however, experiments are technically challenging, time consuming, and expensive to 
perform at the atomistic level. Numerical methods are being used to reduce the 
cost and decrease the time of material and device development. Properly resolving 
these material defects, however, requires full atomistic resolution at the core of the 
defect. Full atomistic resolution throughout a system quickly leads to computational 
intractability for length scales of interest. Multiscale methods like the Quasicontin- 
uum (QC) Method [66, 67, 51, 40, 68], the Bridging Domain Method [78], and others 
were developed to provide atomistic resolution at defects, while coarse-graining far- 
ther away to reduce computational cost. These methods have been used to provide 
information about atomistic defects and structures for a variety of materials. 

In this chapter, we focus on a specific version of the QC method developed in 
[40], which we call the non-local force QC method (NL-QC). This method formulates 
the atomistic problem at the force level and solves a system equations satisfying 
F = 0, where F represents the forces on atoms in the system. The method has 
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currently been developed for simulating atomistic defects in materials with only a 
single species of atoms, for example, pure aluminum, nickel, copper, etc. Some work 
was previously done to extend the method to materials with multiple species of 
atoms (barium titanate and lead titanate) by Kowalewsky in her PhD thesis [42], 
but results were never published and access to the computational implementation is 
not available. Additionally, the work only allowed for simple pair potentials. While 
many body potentials like EAM [16] are conceptually similar to pair potentials, 
their computational implementation can be much more challenging in a multiscale 
framework. EAM potentials, however, provide some of the best models for a variety 
of materials and are essentially required to simulate certain alloys. The research 
presented in this chapter focuses on developing an extension to the NL-QC method 
for multiple species of atoms with no restrictions to interatomic pair potentials. We 


seek to achieve this goal through the following: 


e Developing a mathematical framework for multiple lattices of atoms which 


results in the same formulation for the initial method developed in [40]. 


e Leveraging the existing computational implementation to speed development 


of the multi-lattice implementation. 


e Testing the newly developed method with several patch tests to verify the 


correctness of the code. 


e Analyzing void relaxation in nickel aluminide (B2-NiAl) subjected to volumet- 
ric expansion using EAM potentials [52, 8] to highlight an application of the 


new method. 
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3.3 Methodology 


3.3.1 Overview 


The NL-QC method, like most atomistic to continuum multiscale methods, addresses 
two main computational challenges. The first is the reduction of system degrees of 
freedom. Fully atomistic simulations at length scales of interest in device design can 
quickly approach 107? degrees of freedom or higher. Performing these simulations 
is currently not computationally tractable. The number of degrees of freedom are 
reduced by keeping full atomistics at areas of interest and only tracking a subset of 
atoms elsewhere, called representative atoms in the case of NL-QC. The locations 
of all other atoms in the system are then approximated by shape functions attached 
to the rep atoms, the shape functions used in NL-QC are standard C° continuous 
functions from traditional finite element methods. This coarse-graining immediately 
reduces the number of degrees of freedom drastically and the system size is reduced to 
a computationally tractable level. The second main challenge is the rapid evaluation 
of the force or energy landscape throughout the system. Regardless of the number of 
atoms explicitly tracked with rep atoms, if the force or energy calculation needs to 
be performed at every atom in the system, the computational expense will quickly 
grow. This difficulty is overcome through the use of quadrature points. The use 
of numerical quadrature allows the forces to be evaluated at only a small subset of 
atoms. The NL-QC method uses cluster summation rules located at rep atoms for 
the rapid evaluation of forces. The combination of these two schemes in the NL-QC 
method allows single species systems to be analyzed that could not be handled with 
full atomistic simulations. 

The initial implementation of the NL-QC method took a considerable amount of 
time and effort, especially in regards to the design of data structures and algorithms. 


The code was developed using a combination of fortran, C, and C++ and addi- 
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tionally takes advantage of PThreads for simulation in high performance computer 
environments. In light of the considerable development previously undertaken, we 
seek to utilize a significant portion of the initial single species NL-QC computational 
framework in the extension to multiple species of atoms. We achieve this by treating 
the problem as the union of many single species problems with additional function- 
ality added at a few key locations. We start by reviewing the initial NL-QC method 
and then extending the method to many species of atoms within a properly defined 
mathematical framework. Approaching the problem in this manner allowed the max- 
imum amount of code reuse and drastically reduced the mathematical complexity of 


the formulation. 
3.3.2 Single Species: Problem statement 


We start by considering a crystalline material within a d, dimensional reference 
domain, 2, see Figure 3.1. The material initially consists of periodically repeating 
array of atomistic sites mapped one-to-one with a Bravais lattice, £L. We define a 


Bravais lattice, using basis vectors a;,7 = 1,...,d, as the following: 


L(a) = (X € R?, X =l'a; where ['€ Z',i=1,... _d) ; (3.1) 


Each coordinate, X, in §2 identifies and locates an atom at a unique lattice site. 
We define x as the location of all atoms in the deformed (current) configuration. 
The same bookkeeping system as the reference configuration is used, but all atoms 
are allowed to deform. We initially assume that every atom in the system interacts 
with every other atom in the system through interatomic potentials ranging from 
simple pair potentials like Lennard-Jones to many body potentials like EAM [16]. 
The energy of an atom is now the sum of the interatomic potential with all other 
atoms. The potential energy, F(a), in 2 is just the sum of all atomistic energies 


in the system. We write the total energy of the system, including external loading 
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Figure 3.1: Simple lattice in body 


terms as: 


O(a) = E(x) + 6°"(2). (3.2) 


@*'(x) includes all contributions from Neumann or Dirichlet boundary condi- 
tions, body loads, and/or any other external loading device. We seek to find local 


minimizers of (a), which correspond to equilibrium configurations. 


min &(a) (3.3) 


x 


3.3.3 Single Species: Reduction of system degrees of freedom 


In fully atomistic systems, the number of degrees of freedom quickly reaches compu- 
tationally intractable levels for even modest sized domains. Modeling the complete 
atomistic system is not possible at this time. Thus, we reduce the number of degrees 


of freedom in the system by selecting a subset of atomistic sites 1, C l with atomistic 
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locations X;, C X and a, C 2 to explicitly model. We call this subset of atoms 


‘representative atoms’. We now seek a minimization over the representative atoms. 


min &(a2) (3.4) 


Lh 


(x) is still dependent on the position of all the atoms in (2, but we are only explicitly 
tracking the representative atoms. We introduce a triangulation, it does not need 
to be structured, 7, over the representative atoms X;, see Figure 3.2, where the 
dark blue atoms are representative atoms. The top zoom out is a coarse-grained 


mesh, while the bottom zoom out is a fully resolved mesh. Shape (interpolation) 
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Figure 3.2: Simple lattice with full atomistic and coarse-grained resolution 


functions defined on the simplices of the triangulation are then used to approximate 


the locations of all non representative atoms in the system. 


v= 5), (.X(1),2%) (3.5) 
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The above equation maps atomistic sites to their location in the current configuration 
through the, S;, shape functions. In principle, any form or type of shape functions 
can be used, though throughout this chapter we define S, to be standard piecewise 
linear shape functions. Given this definition of S;, the location of all atoms in 2 is 


only dependent on a). &(-) can now be written as a function of ap. 
3.3.4 Single Species: Reduced equilibrium equations 


The full energy minimization problem from Equation 3.3 can be rewritten in terms 


of the reduced energy: 


min &(x,),). (3.6) 

Lp, 
The reduced equilibrium equations, which correspond to local energy minimizers are: 
Fil @y,) = Oe Fe), 24) = 0. (3.7) 


F;,(a,) is the total force at representative nodes a,. x = S;,(X(l),x,) maps the 
atomistic locations in the reference to the current through the representative atom 
positions and triangulation shape functions. S,(F(a),x;,) maps the forces on a to 
forces on representative atoms a, through the triangulation shape functions. The 
forces at a given atom a; are: 

F (te) =O g(t) = S- W,(r)?.2,, where (3.8) 

@pEP; (xi) 

W(r) is the known interatomic potential for a given material and r is the distance 


between atoms 7 and j, |a; — a;|. Additionally, P, is defined as follows: 
PAD) = {tae |e — Oi| = Ay t F797) whete (3.9) 


R, is the potential cutoff radius. 
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Essentially, the forces are calculated at all atoms in (2 and distributed to rep- 
resentative atoms, 2,, via the same linear shape functions used to determine the 
location of the non representative atoms in (2. The reduced equilibrium equations 
currently still require calculating the force at every atom x in (2, where the force 
is calculated by the full lattice sum of an atom interacting with every other atom 


within a cutoff radius. There are two difficulties that need to be overcome: 


1. Evaluating the forces at every atom is prohibitively expensive. Quadrature 
rules need to be introduced to evaluate the forces at selected quadrature points 
to approximate the full system. We use the same cluster summation rules 
centered at representative atoms from the initial NL-QC development [40]. 


The specific form of the cluster summation rule will be explained below. 


2. In the full atomistic calculation every atom interacts with every other atom 
in the system. We decrease the number of calculations by assuming that all 
interatomic potentials are short-ranged and can be approximated with an ef- 
fective cutoff radius, R,. Any atoms outside the cutoff radius will not effect 
the evaluation of the forces. When we say short-range, we mean that the in- 
teratomic potential has to decay at least as fast as r~*, where r is the distance 
between two atoms. This assumption eliminates the modeling of Coulombic 
and Coulombic-like interactions and restricts the class of materials that can be 
simulated. Previous work has extended some QC methods to handle Coulom- 
bic interactions [49], while work is ongoing to extend the NL-QC to Coulombic 
interactions, thereby relaxing the short-ranged interatomic potential assump- 


tion. 
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3.3.5 Single Species: Summation rules 


As stated above the full force calculation, essentially \7¢9 F(x), is very compu- 
tationally expensive and needs to be reduced for realistic length scales of (2. This 
reduction is achieved through the use of quadrature. Forces are calculated at a sub- 
set of atoms, x, C {a} and may or may not have any relation to the representative 
atoms a,. These forces at quadrature points are then multiplied by a weight that 
is essentially the number of atoms that the quadrature point approximates and is 
analagous to volume in a classical continuum mechanics setting. The reduced equi- 
librium equation can now be written as: 
Flite is Ss Q(&q)Sn(F (x,),£n) = 0, where (3.10) 
LqgC LED 
Q(a,) is the quadrature weight at 2,. The above equation is the general form 
of quadrature, but many different types of quadrature rules can be selected. For 
example, the natural first choice for quadrature points would be x, = x;,, where the 
quadrature points are the exact same as the representative atoms. Knap and Ortiz 
[40] showed that this choice of summation rule suffered from rank-deficiency. They 
overcame this difficulty by using a cluster summation rule centered at representative 
atoms. The cluster summation rule evaluates forces at atoms within a given cluster 
radius of a representative atom and uses those forces to approximate the total force 
at a representative atom. The reduced equilibrium equation for this specific rule 
looks like: 
Fr(an)* S> Q(an) | S > Sh(F(a-),an)| = 0. (3.11) 


ap,cxEeQ @-EC(xp) 
C(ap) is a cluster of atoms defined by C(x;,) = {a : |w — a),| < R-}, where R, is the 
cluster radius centered at a,. F(a.) is calculated by summing the force over all of 
the atoms that a given cluster atom 2, interacts with. In a triangulation for a real 
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problem, there will be areas of high and low representative atom density. In areas of 
low density, which correspond to a coarse-grained approximation of the atomistics, 
clusters will not overlap and Equation 3.11 is used as expected with no issues. Con- 
versely, in areas of high density, specifically where triangulation is reduced to every 
atom being a representative atom, the cluster summation rule evolves directly to 
a full force calculation. This situation corresponds to a fully resolved, completely 
atomistic area. The intermediate area between a coarse-grained and a fully resolved 
triangulation can have overlaps between clusters. The question in this case is which 
cluster should an atom a belong to? This question is resolved by choosing the nearest 
representative atom, in the case that two or more representative atoms are the same 
distance, the conflict is resolved randomly. Figure 3.3 highlights the three different 


scenarios. 


Figure 3.3: Simple lattice with multiple cluster scenarios. Top Zoom: Coarse- 
grained no cluster overlap. Middle Zoom: Coarse-grained with cluster overlap that 
must be resolved. Bottom Zoom: Fully resolved cluster is a single representative 
atom. 
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3.3.6 Single Species: Adaptive meshing 


In the previous section, we discussed various mesh refinement levels, but never dis- 
cussed how those were obtained. At the beginning of a simulation, a generic mesh is 
known a priori, tailored to the specifics of the problem. A fully resolved mesh is spec- 
ified around atomistic defects and is slowly coarse-grained with increasing distance 
from the atomistic defect or defects. It is unknown if an initial mesh has enough 
atomistic resolution for a given applied loading, so an adaptive meshing algorithm is 
used. This algorithm refines the mesh if the strain is above a given tolerance level. 
In practice, the tolerance level is set so that crystallographic slips of a full Burgers 
vector are resolved with full atomistic resolution. Specific details of the adaptive 


meshing algorithm and tolerance parameter can be found in [40]. 
3.3.7 Multiple species: Problem statement 


In the previous sections, we outlined the single species NL-QC method. We now seek 
to extend this method to multiple species of atoms by treating the problem as the 
union of multiple single species problems. We again start by considering a crystalline 
material within 2, except now there are multiple species of atoms, N*, denoted by 
Greek letters (a, @, etc.). The same definition of a simple Bravais lattice is used, 


except a few modifications are made. 


L°(a*) = (X* ERY, X% = [a2 + b® where [** € Z',i=1,...,d). (3.12) 


1. We introduce a shift vector, b®, that denotes the relative displacement between 
lattices. This shift vector establishes a global coordinate system. Essentially, 
each lattice has its own local coordinate system determined by I°%’a%, that is 
transformed to the global coordinate system through the shift vector. Typi- 
cally, we take the 0" lattice to have the same local coordinate system as the 
global coordinate system by setting b° = 0, though this choice does not have 
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to be made. 


2. We allow each simple lattice to have different lattice vectors | and lattice con- 
stants a. These different vectors and constants are very important for modeling 
different material phases. A good example of this occurs in NiAl, which can 


form in many different phases including B2-NiAl and L1,2-Ni3Al. 


B2-NiAl] is perhaps the simplest case of a multi-lattice and can be thought of as 
a cubic lattice of Ni and a cubic lattice of Al offset by half of the lattice constant 
a in ad dimensional space. In this case, the lattice vectors and constants are the 
same, 1° = I' and a° = a!, where 0 and 1 denote Ni and Al. The shift between the 
two lattices can then be set by defining b? = 0 and bj = I'"a}/2. 

The L1,-Ni3Al phase is slightly more complicated. This phase can be thought of 
as a cubic lattice of Al intersected with an fcc lattice of Ni. In this case the shift 
vectors b® are important, but additionally, the lattice vectors /* and lattice constants 
a* are different and describe a cubic and fcc lattice respectively. The different lattice 
vectors and constants allow much more complicated phases to be easily modeled and 
are pivotal for being able to simulate many different phases. 

We now define a complex Bravais lattice as L, where all of the atomistic locations 


in the reference coordinate system are denoted by X. 
L=2UL Uesul,. X= Xo UX UX (33) 
For illustration purposes only, N° = 2 is chosen and shown in Figure 3.4. The location 


of all atoms in the current configuration, energy in 2, and energy minimization 


equations proceed as expected: 


X=2Uaelu---Ue™, (3.14) 
O(X) = E(X) + 6" (2X). (3.15) 
min ®(¥) (3.16) 


Offset by shift 


Figure 3.4: Multiple Species Lattice 


3.3.8 Multiple Species: Reduction of system degrees of freedom 


We now seek to reduce the number of degrees of freedom in {2 by selecting a subset 
of each a lattice to be representative atoms. These representative atoms will only be 
used to interpolate and simulate atoms within the same lattice, i.e. Ni representative 
atoms would not be used to interpolate Al atoms in NiAl. Thus, we select a subset of 
atomistic sites [’ C 1° with atomistic locations X;* C X° and a? C 2x to explicitly 
model. The full set of representative atoms is now V%, = af Ua, U---Uah”. 
Triangulations, 7,°, are introduced for each lattice £° using representative atoms 
X;, see Figure 3.5. A set of piecewise linear shape functions are defined over each 


triangulation to interpolate positions of non representative atoms of the same species. 
a = Sy (XS (1), ey) (3.17) 
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There are several important points to note: 


e Triangulations need not be structured. 


e Triangulations for different species of atoms do not need to have the same 


structure. 


e There is no limit on the number of species in a material, and thus no limit on 


the number of triangulations. 


e Because of the separate triangulations, adaptive meshing stays the exact same 


as in the single species case for each lattice. 


Offset by shift 


Figure 3.5: Multiple Species Lattice with Multiple Triangulations 
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With the equations to locate every atomic location in {2 now defined, the minimiza- 


tion of the energy occurs over all of the representative atoms. 
min &(X,) (3.18) 
Xn 


3.3.9 Multiple Species: Summation rules 


We now seek to solve the reduced equilibrium equations for (2, this amounts to 


simultaneously solving the equilibrium equations for all N° species of atoms. 
F(a) =0, Fy (vj) =0, ..., Fi” (az) =0 (3.19) 


We use the same cluster summation rules for a single species of atom for each a‘” 


species of atom to reduce the force calculation. 


Fran) DS) Qa) | Dd) Si(Fe(we), wf) | = 0. (3.20) 
a? Care ere (ae) 

Until this point, all calculations for the a” species have been treated completely 
separately. That distinction changes with the calculation of the forces on an atom. 
An individual atom ’sees’ other atoms through interatomic potentials, while the po- 
tentials may change depending on which species are interacting, fundamentally there 
is no unique distinction between atoms of different species. This is the main differ- 
ence between multi-lattice NL-QC and simple lattice NL-QC, during the calculation 
of forces, atom a sees all other atoms within a cutoff radius, regardless of species. 


Thus, the force calculation is now: 
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W°"(r) is the interatomic potential between a and 6 species of atoms and r is the 


distance between atoms 7 and 7, xi —ax°%|. Additionally, P® is defined as follows: 
p= {ae Jai = 07 |<S ht y i o= 8}, where (3.22) 


R, is the potential cutoff radius centered at x*. Figure 3.6 highlights the formation 


of clusters and how atoms interact in the force calculation. 


Figure 3.6: Multiple Species Lattice Clusters Top Zoom: Light blue and green 
atoms are in a cluster of dark blue and green representative atoms respectively. 
Bottom Zoom: The calculation of forces at dark blue and green atoms sums over 
interaction with all atoms within a cutoff radius regardless of species. 


Notes 


e The notation for multiple species of atoms evolves exactly to the existing for- 
mulation for the single species NL-QC, if the number of species is taken as 1. 
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This evolution further ensures the consistency of the multiple species method 


with the initial formulation 


e An initial error estimate of the single species cluster summation rule was con- 
ducted in [40] based on different sizes of cluster radii, calculation of cluster 
weights, and other factors. A subsequent analysis conducted in [46] showed 
that the cluster summation rule introduced error into the model. Additional 
quadrature rules were proposed for the NL-QC method to reduce this error 
[30, 34], however, these rules have their own issues. At this time, it is our opin- 
ion that there is no clear ’best’ quadrature rule currently available in literature. 
Research is ongoing in the QC field to develop better quadrature rules. For 
these reasons, we chose to extend NL-QC to multiple species using the same 
cluster summation rule developed in [40]. If a ’best’ quadrature rule becomes 
available, conversion of the code and details of multiple species NL-QC to this 
quadrature rule can be undertaken following the same procedure illustrated 


above. 


3.4 Computational implementation 


3.4.1 Overview of single species NL-QC code 


The single species NL-QC code was initially developed using a combination of For- 
tran, C, and C++. Additionally, extensive work was done on implementing advanced 
data structures, enhancing calculation speeds using local processor caches, and par- 
allelizing algorithms using PThreads. The code starts by initializing data structures, 
building neighbor lists, computing forces on nodes, and then minimizing the forces 


using a nonlinear conjugate gradient algorithm. The general flow of the code is 
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below: 
InitializeDataStructures(...); 
BuildNeighborList(...); 
ComputeForces(...); 
while(not converged) { 
RebuildNeighborListIfNecessary(...); 
ComputeForces(...); 


UpdateRepAtomLocationsWithCG(...); 


} 


OutputResults(...); 
3.4.2 C++ class extension of Quasicontinuum to Quasicontinua 


In developing the extension of NL-QC to multiple species, our goals was to reuse as 
much of the existing code as possible. We achieved this by creating a C++ class for 
a single Quasicontinuum, which holds all of the data structures for a single species 
or lattice of atoms. Essentially all of the data structures for the original code are 
replicated and initialized each time we instantiate a new Quasicontinuum lattice. For 
example, if we wanted to simulate NiAl, we would instantiate a Ni Quasicontinuum 
lattice and an Al Quasicontinuum lattice, where each Quasicontinuum lattice keeps 
track of its own data. A Quasicontinua class was then created to hold all N* Qua- 
sicontinuum instances and enable easy access to each Quasicontinuum. Additional 
C++ classes were implemented for communicating necessary data between different 
lattices. 

One of these additional C++ classes is the BuildCrossNeighborList class, which 
calculates and stores all atomistic neighbors of a site, regardless of which atomic 


species or lattice the neighbor belongs to. The BuildCrossNeighborList class stores 
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the neighbor sites and locations in separate hash tables for quick searches and min- 
imal duplication of data. A second class was created to allow implementation of 
EAM potentials. EAM potentials are different from traditional pair potentials, like 
Lennard-Jones, in that they are a many body potential. Equations 3.23 and 3.24 


show the form of the energy in a system of atoms. 


1 _ 
Brn = 5 S> Yass (Tis) + > _ Fa, (Ai), where (3.23) 
1,5 (J#t) a 
i and j are individual atoms, a; is the species of the 7” atom, a; is the species of the 


j atom, Fy, is the embedding energy, and is the total charge density of an atom. 


pi = Seg (riz), where (3.24) 
iAj 

Pa; is a known function of charge density. The many body nature of the potential 
comes from the fact that to calculate the embedding energy, the charge density at 
an atom must be known. However, the charge density of an atom is dependent on 
all other atoms within a given cutoff radius all of which must be known before the 
calculation of forces. We created the EAM class to precompute the charge density, 
based on CrossNeighborList data, at all necessary atoms before the density is needed 
to compute forces and energies. Finally, a ComputeCrossForces class was created, 
which takes neighbor list and EAM data as an input, and calculates all forces on 
atoms in clusters and subsequently forces on representative atoms. All of these classes 
were developed using PThreads for seamless integration into the existing code for 


high performance computing environments. The computational flow of the extension 
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of NL-QC to multiple species of atoms is below. 

for(i=0;i< N*;++i){ 
iQuasicontinuum = Quasicontinuum::InitializeDataStructures(...); 
Quasicontinua::insert(iQuasicontinuum); 

} 

for (i=0;i< N*3++i){ 
iQuasicontinuum = Quasicontinua::getQuasicontinuum(i); 
iQuasicontinuum.BuildNeighborList(...); 
iQuasicontinuum.BuildCrossNeighborList(...); 
iQuasicontinuum.EAM(...); 

} 

for (i=0;i< N*3++i){ 
iQuasicontinuum = Quasicontinua::getQuasicontinuum(i); 
iQuasicontinuum.ComputeForces(...); 


iQuasicontinuum.ComputeCrossForces(...); 
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while(not converged) { 

for(i=0;i< N*;++i){ 
iQuasicontinuum = Quasicontinua::getQuasicontinuum(i); 
iQuasicontinuum.RebuildNeighborListIfNecessary/(...); 
iQuasicontinuum. RebuildCrossNeighborListIf{Necessary(...); 
iQuasicontinuum.EAM(...); 

} 

for(i=0;i< N*;++i){ 
iQuasicontinuum = Quasicontinua::getQuasicontinuum(i); 
iQuasicontinuum.ComputeForces(...); 
iQuasicontinuum.ComputeCrossForces(...); 

} 

for (i=0;i< N*;++i){ 
iQuasicontinuum = Quasicontinua::getQuasicontinuum (i); 


iQuasicontinuum. UpdateRepAtomLocationsWithCG (...); 


} 
for(i=0;i< N*;++i){ 
iQuasicontinuum = Quasicontinua::getQuasicontinuum(i); 


iQuasicontinuum.OutputResults(...); 
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3.4.3 Verification of implementation 


We verified the new code with several simple patch tests. The first test, verified the 
correct calculation of 0 forces and charge densities in undeformed B2-NiAl using the 
potential developed in [52]. In the next test we applied a simple affine shear deforma- 
tion to the system. The results again produced 0 forces in the system, as required. 
Finally, we applied the same affine shear deformation, but this time only to the 
boundary nodes. The system was then allowed to minimize using the implemented 
nonlinear conjugate gradient algorithm. The final equilibrium atomistic positions 
corresponded exactly to the expected affine deformation, to within numerical toler- 
ance. The passing of these simple tests gives us a confidence that the implementation 


is correct. 
3.5 Numerical simulation of deformation around void defects in NiAl 


We use the method to simulate the deformation around various sized voids in NiAl 
subject to volumetric expansion. NiAl is modeled in its B2 phase with an EAM 
potential [52, 8]. We show a typical mesh used in Figure 3.7. Near the void, full 
atomistic resolution is used with coarse-graining occurring farther away from the 
defect. Void sizes from 2-10 Angstroms are simulated subject to load increments of 
volumetric expansion. We apply the deformation to the system, let the atoms relax 
and then apply the next loading. Figure 3.8 and Figure 3.9 show the centrosymme- 
try parameter around a 2 Angstrom void (in this case only a vacancy) after initial 
relaxation and at a volumetric expansion of 24 percent. Figure 3.10 and Figure 3.11 
show the centrosymmetry parameter around a 4 Angstrom void after initial relax- 
ation and at a volumetric expansion of 16.5 percent. Figure 3.12 and Figure 3.13 
show the centrosymmetry parameter around a 6 Angstrom void after initial relax- 


ation and at a volumetric expansion of 15 percent. Figure 3.14 and Figure 3.15 show 
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Figure 3.7: Typical mesh with coarse-graining far away from void with full atomistic 
resolution 


the centrosymmetry parameter around a 8 Angstrom void after initial relaxation 
and at a volumetric expansion of 13.5 percent. Figure 3.16 and Figure 3.17 show 
the centrosymmetry parameter around a 10 Angstrom void at after initial relaxation 


and at a volumetric expansion of 12 percent. 


81 


Pseudocolor 
Var: csym 
5.953 


4.465 
2.977 
1.488 


0,000 
Max: 5.953 
Min: 0.000 


Figure 3.8: Centrosymmetry parameter around vacancy after initial relaxation. 
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Figure 3.9: Centrosymmetry parameter around vacancy at volumetric expansion of 
24 percent. 
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Figure 3.10: Centrosymmetry parameter around a void of 4 Ang. after initial relax- 
ation. 
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Figure 3.11: Centrosymmetry parameter around a void of 4 Ang. after volumetric 
expansion of 16.5 percent. 
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Figure 3.12: Centrosymmetry parameter around a void of 6 Ang. after initial relax- 
ation. 
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Figure 3.13: Centrosymmetry parameter around a void of 6 Ang. after volumetric 
expansion of 15 percent. 
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Figure 3.14: Centrosymmetry parameter around a void of 8 Ang. after initial relax- 
ation. 


Pseudocolor 
Var: csym 
16. 


12,37 


8.249 


4.125 


0.000 
Max: 16.50 
Min: 0.000 


Figure 3.15: Centrosymmetry parameter around a void of 8 Ang. after volumetric 
expansion of 13.5 percent. 
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Figure 3.16: Centrosymmetry parameter around a void of 10 Ang. after initial 
relaxation. 


Pseudocolor 
Var: csym 
25.49 


19.11 


12.74 


6.371 


7.414e-06 
Max: 25.49 
Min: 7.414¢-06 


Figure 3.17: Centrosymmetry parameter around a void of 10 Ang. after volumetric 
expansion of 12 percent. 
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3.6 Discussion 


In this chapter, we extended the non-local force QC method to multiple lattices. We 
treated the problem as the union of many simple lattice problems and then allowed 
the force calculation to communicate relevant information between the lattices. The 
method was verified with a variety of simple tests after implementation in a high 
performance computing framework. We used the method to simulate void relaxation 
in B2 NiAl for multiple void sizes from a vacancy to lnm. We intend to continue 
loading these voids until dislocation nucleation occurs from the void to understand 


how void size effects plasticity. 
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A 


Non-local quasicontinuum method with long-range, 
real-space electrostatics 


4.1 Abstract 


We develop a method for simulating ionic solids with long-range electrostatic in- 
teractions, without assuming periodicity. The new methodology can be described 
as a combination of the Reaction Field and Fast Multipole methods. The method- 
ology we use allows for the evolution of a body potentially with defects (domain 
walls, vacancies, etc.) subject to electrostatic and mechanical loadings. Additionally 
we allow explicit handling of multiple types of electrical boundary conditions. The 
problem we are interested in is essentially calculating the electric field at a given 
location from a multi-lattice of charges. Generally speaking, we want to break up 
the electric field calculation into two components. The first component is an exact 
summation region near the location of interest. The second component includes all 
parts of a body outside of this region that are away from defects. The problem is 
too computationally expensive to perform exact summations everywhere, so we cal- 


culate an approximated electric field from the far field region using ideas from the 
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Fast Multipole Method. 

In this chapter, we derive the method by performing a Taylor series expansion 
of the standard Coulombic kernel for charges. We start by considering a single unit 
cell of charges, where the net charge of the unit cell is zero. We work only with 
the dipole expansion term and consider all higher order moments as bounds on the 
error of the expansion. A summation over a large number of unit cells in a body is 
performed. The dipole term is transformed into a standard divergence term, while 
the error term is integrated over shells of expanding unit cells centered at the point of 
interest. An error bound for the calculation of the electric field at a point of interest 
is analytically derived and numerically modeled for GaN, though the analysis can be 
extended to any ionic solid material. We then use the error bound as a criterion to 
determine the size of the exact summation region around an evaluation point for the 
electric field. Our work has been implemented within an existing Quasicontinuum 
Method code for multiscale modeling and high performance computing. 

Keywords: Quasicontinuum (QC) method; electrostatics; Fast Multipole Method(FMM); 
Reaction Field(RF) method; ionic solids 


4.2 Introduction 


Electrostatic (Coulombic) interactions play a prominent role in the structure of all 
materials. These interactions fundamentally take place at the quantum length scale 
and occur between electrons and protons. Modeling materials at this fundamen- 
tal length scale is quite challenging, even with sophisticated methods like density- 
functional theory (DFT). DFT simulations are typically comprised of hundreds of 
atoms, but occasionally are used to model systems with thousands or even tens of 
thousands of atoms, though these simulations are rare. These simulations struggle to 
provide detail at the larger length scales important in engineering, where the num- 
ber of atoms in a system can be orders of magnitude higher than DFT can handle. 
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Conversely, classical continuum methods provide coarse-grained information at the 
engineering design length scale, but cannot fully resolve the underlying structure of 
materials. Atomistic methods like molecular dynamics [23] and multiscale methods 
like QC [68, 51] are in between the two in terms of both the largest accessible length 
scale and highest resolution of the underlying material structure. The material in- 
put parameters in atomistic methods can be thought of as the interatomic potentials. 
Interatomic potentials essentially approximate the underlying energy and force inter- 
actions between atoms ranging from bonds between pairs of atoms (Lennard-Jones, 
Buckingham, etc. [23]) to many body interactions (EAM[16, 17], MEAM[5], etc.). 
Many interatomic potentials do not explicitly include Coulombic terms, but instead 
include them in an approximate way. The underlying structure of the material is 
essential to this choice. 

The Coulombic interaction kernel is shown in Equation 4.1, where q; and q; are 
the effective charge of two different atoms and r;; is the separation between the two 


charges. 


Gd; 
= 4.1 
Aneori; ( ) 


This interaction energy is long-ranged, because of its slow 1/r decay and can be a 
divergent series when summed over many interacting atoms depending on the under- 
lying structure. In ionic solids like (GaN, Ceria, BaTiOs3), the situation is better and 
the Coulombic energy effectively decays as 1/r?, because of the orientation of atomic 
charges in the crystal structure. Interatomic potentials for ionic solids explicitly in- 
clude these Coulombic interactions. Calculating the energy in a system then requires 
adding the interactions between all atoms in the system, resulting in a conditionally 
convergent sum due to the effective 1/r? decay. Metallic solids, however, are differ- 
ent because the Coulombic energy decays as 1/r° due to the different orientation of 


atomic charges in the crystal structure. In this case, most interatomic potentials do 
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not explicitly include Coulombic interactions. Instead the Coulombic energy is in- 
cluded with the other short-range interactions. This is valid because the Coulombic 
sum is absolutely convergent in metals and a cutoff radius typical for short-range po- 
tentials can be applied. A more detailed description of why electrostatic interactions 
are long-ranged can be found in [49]. 

In this chapter, we focus specifically on ionic solids, where the electrostatic inter- 
action is explicitly included in an interatomic potential. While the problem starting 
from the quantum length scale has already been simplified with an effective charge 
on atoms instead of handling individual protons and electrons, performing the full 
summation is still quite challenging. The full atom-atom interaction summation is 
not feasible due to its large computational cost. In its place, many methods have 
been developed to approximate the summation. One class of methods are the Ewald 
[21, 23] and Ewald-like summations [22, 23]. These methods break up the full sum- 
mation into a real-space component and an imaginary space component, both of 
which are fast to calculate. The real space component adds the Coulombic interac- 
tions from atoms within a small cutoff radius r, while the imaginary space component 
adds the interactions from all atoms outside the cutoff radius. These methods are 


quick and effective, but they suffer from a few drawbacks. 


e Boundary conditions beyond simple (cube, sphere) geometries are difficult to 


handle. 


e The imaginary space component assumes some form of periodicity to perform 
the summation. If the far field atoms are varying in position from the assumed 


periodicity, then this assumption does not capture that variation. 


Another prominent method for electrostatics is the fast multipole method (FMM) 
[28, 7, 14]. This method performs a Taylor series expansion of the Coulombic kernel 
about a cluster of charges. The method then calculates the higher order moments 
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from the cluster and approximates the net interaction with the resulting moments. 
The error in the method can be controlled based on the distance of a cluster from 
the point of interest and the number of moments included in the approximation. 
FMM, however, is not as widely used in the atomistic modeling world as other fields. 
Some of the challenges are its complex algorithms and subsequent high degree of 
computational implementation difficulty, especially in a high performance computing 
environment. 

The last prominent method that we want to highlight is the Reaction Field (RF) 
method [73, 27, 6]. This method is prominently used in computational chemistry for 
molecules surrounded by a homogeneous liquid or solid. The Coulombic interactions 
in this method are included exactly within a cutoff radius around a molecule or 
molecules. Outside of this cutoff radius, a far field term arises by assuming an 
electrically active homogeneous medium typically a liquid or a solid. This medium 
then acts on the molecules within the cutoff radius via a constant electric field arising 
from continuum level calculations. The exact values of the field are calculated by 
solving Maxwell’s electrostatic equation in the body containing the homogeneous 
medium. The far field domain has two important boundaries, one on the exterior 
and one on the interior around the exact summation region that result in the applied 
field. The interior boundary is typically taken as a sphere, though it does not need 
to be. 

In the atomistic mechanics field, the Ewald and Ewald-like methods are most 
prevalently used, but introduce the serious limitations stated above. We seek to 
develop a method that does not have these limitations, which we achieve in a new 
method that can be described as a combination of FMM and the RF method. Our 
work in this chapter extends our initial research ({49], Chapter 2) on modeling mate- 
rials with Coulombic interactions. The enhanced method developed in this chapter 
can simulate 3-D problems, is fully non-local, unlike our previous work, and is im- 
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plemented in a high performance computing framework. 


4.3 Method 


4.3.1 Preliminaries 


We start by assuming an atomistic body, 92, see Figure 4.1. (2 is a perfect single 
crystal with multiple species of atoms per unit cell. The location of unit cells are 
denoted by X in the reference configuration and x in the current configuration, while 
the locations of specific atoms within a unit cell are denoted by y*. Each unit cell is 
referenced by a unique triple integer Z € Z® that locates the specific unit cell within 
a lattice. We assume an effective charge on each atom that is denoted by q*, where 
a is the species of atom in the unit cell. The charge on an atom does not need to 
be discrete, but we assume that it is and extend to a continuum charge density later 
on with a Dirac delta function centered at all y® atomic locations. Additionally, all 


unit cells of atoms must be charge neutral and have a length scale of I. 


Sig? =0 (4.2) 


Figure 4.1: Body (2 with atomistic charges and charge neutral unit cells. 
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4.3.2 Quasicontinuum Method 


The QC method is a powerful tool for solving large atomistic systems. The method 
was originally developed by Tadmor, et. al. in 1996 [66]. Since then multiple versions 
of the method have been developed, more information on the specific methods can 
be found in these review articles [68, 51]. In our previous work with electrostatic 
interactions, we used a purely Cauchy-Born variant of QC [49]. This version of QC 
is valid in large elements, where the deformation gradient is slowly varying. The 
version is not valid near defects or areas of high deformation, which we are mainly 
interested in due to the interesting physics occurring in these regions. For this work 
we use the non-local cluster force QC [40], specifically a variant for multi-lattices 
[42, 50] that is valid near areas of high deformation. 

We will now give a brief overview of the QC method, for more in depth details, 
see [40, 42, 50]. In this version of QC, we are interested in finding local minimizers of 
the energy, (a + y®), in a body Q, see Equation 4.3. These minimizers correspond 
to equilibrium configurations of (2 subject to any defects or external loadings in the 


system. In the remainder of this section, we will use the notation #° = 2+ y’°. 


min @(a") (4.3) 


As an atomistic system becomes large, explicitly tracking all atoms in (2 is not com- 
putationally possible, so we introduce a coarse-graining approximation in slowly de- 
forming regions. The coarse-graining is achieved by only tracking a subset of atoms, 
which we call representative atoms, ah. In areas with defects or external loadings, 
every atom 2° is a representative atom a, to achieve full atomistic resolution. The 
regions of high deformation are initially chosen at the beginning of a calculation, but 
adaptive meshing has been implemented and can be used to adjust the fully resolved 
regions, where necessary. In regions far away from the high deformation areas, we 
only select a small subset of atoms to be representative atoms. The locations of all 
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other atoms are then interpolated using standard shape functions, S? , based on the 


representative atoms. 


2? = So (x*(2), w}) (4.4) 


In our multi-lattice case, each atom within a unit cell will be interpolated from a set 
of shape functions, S? , and representative atoms a’ of the same § species. Hence, 
there are mazx(3) sets of representative atoms and triangulated meshes that shape 
functions, S? , are built from. 

We can now write the equilibrium configurations of (2 in terms of the represen- 
tative atoms. 


min P(x) (4.5) 


rp, 


Local minimizing states of these equilibrium configurations correspond to the forces, 
F*(#), equaling zero in 2. However, the computational expense is again too high 
to solve for the forces for all atoms in a system. Thus, the forces are only explicitly 
calculated on representative atoms, where an approximated force is calculated from 


a cluster summation rule [40]. 


Fale + O@ SF aa) (4.6) 


aw EC(xy) 


C(x?) denotes a cluster of «° atoms centered around a representative atom, a’. 
Q(x") is the summation weight on the cluster atom. Equation 4.6 essentially calcu- 
lates a weighted force on a representative atom from a cluster of atoms, this effective 
force can then be used to find equilibrium configurations in 2. Each atom can only 
belong to a single cluster to meet the summation rule requirements in [40] and only 
atoms of the same ( species are used in a cluster around a’. Thus, in areas of 
high deformation, where every atom is a representative atom, the force calculation is 
exact and essentially devolves into molecular statics. Away from these regions, the 
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force is only approximated. We now use a conjugate gradient algorithm to solve for 


zero forces at all representative atoms. 
Fi (x;) = 0 (4.7) 


The only thing still needed to perform the minimization is a way to calculate 
forces on a given atom, F°(a"). We start by breaking up the F°(a*) calculation into 


two components, the typical short-range interactions and the long-range interactions. 


F(a") = Forore(@) + Fe (") (4.8) 


s electrostatic 


QC methods typically calculate F°(x*) using interatomic potentials that are only 
short-ranged with a cutoff radius. Long-range interactions have been included in 
non-local QC methods, but to the best of the authors’ knowledge only with Ewald 
sums [42]. We will not provide any more explanation for calculating short-range 
forces from interatomic potentials, because they have been extensively covered in 
the mechanics literature. Instead we focus specifically on including the long-range 


electrostatic interactions, which we develop in the following section. 
4.3.3 Electrostatic Summation Overview 


Calculating the electrostatic force on an atom requires knowing the electric field E 
and multiplying it by the charge on the atom q’. In this section of the chapter, we 
occasionally use indicial notation where a; denotes a vector and 7 runs over | to the 


number of spatial dimensions of 2. 


ieee Cl) = q’ E,(x*) (4.9) 
We use a standard definition of the electric field, which is related to the electric 
potential by F = —V@¢. The force on a 8 atom can now be written as the explicit 
sum over all a atoms € 2 interacting through the standard kernel. 

B 


g(a, — % — yt) 
Dot) = qr by >, (4.10) 


F' 
Areo|x® — a! — y*|8 


From here, we will focus on calculating the electric field at any given point, x. 
4.3.4 Contribution from Single Unit Cell 


We start with the contributions from a single unit cell as seen in Figure 4.2. The field 


@® l 


Figure 4.2: Schematic of unit cell geometry contributing to the electric field at a 
given point x”. 


can be written in terms of an integral instead of a sum, assuming that each discrete 


charge is approximated by a dirac delta function. The resulting charge distribution 


is p(y). 


= | ply?) (ae _ xi ZZ ys) dy (4.11) 


We are interested in the contribution to the electric field with respect to two com- 
peting length scales |a° —a’|, the distance between the evaluation point and the unit 


cell, and / the characteristic length of a unit cell. We choose the following scalings 


to evaluate these competing length scales, with LJ representing a cell as a unit cube. 


y; = ly;, hence d®’y =I d (4.12) 


p(y) = p(y)/t (4.13) 


The electric field can now be written in terms of the scaling parameter J. 


cel l?p y A a x; = lye 
Are Ee" (ar) = ee aE dg (4.14) 


We now perform a Taylor series expansion of the kernel function about y = 0, 
essentially the multipole expansion commonly used in FMM. We only show the first 
3 terms, but each higher order term adds an additional scaling of / and an additional 


derivative on the kernel. 


Boe 3 (af - 2) (xf - 21) er sal) 
ce ~/ ty Th ~ d g ‘ ‘ 
Aneg ES "(g) y : ?p(y) (a =5 a'|3 yl |x? = an! |® 
1 eee 1 
Ee 41 
+ 5h;G0l? Vi jVe (= — =) + ) d (4.15) 


Note that each of the 3 terms are simply the n™ derivatives, V"(1/r), of the kernel 
function, where r = |x? — a’|. The first term immediately goes away because we 


assume that the unit cell is charge neutral. 


[ a@aa=o (4.16) 


We now use the polarization density in its customary form to modify the first term. 


fd 


All factors of J cancel in the scaled version of the polarization density. Thus, the 
unscaled and scaled polarization are identical. This is the reason why the scaling for 


the charge density was initially chosen as p(y) = p(y)/l. The electric field now can 
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be written in terms of the polarization density. 


: ( 
cell (_,B\ 3 
dren Ee" (w*) ~ p,(") aa 


1 1 
+f 5 Vit! ViV Ve ( =) dd+... (4.18) 


ja? — 


All terms beyond the dipole correspond to higher order multipoles that could be 
explicitly included in the calculation. The FMM does exactly that and controls the 
error in the electric field calculation by including higher order multipole terms as 
necessary. The number of terms to be included depends on the distance between 
the source cluster and the evaluation point namely |a° — a’|. We choose to only 
include the dipole term in the calculation of the electric field and consider all higher 
order terms as a bound on the error in the electric field calculation. As the distance 
|2° — x’| decreases, the error in the electric field has the potential to grow and make 
this approximated calculation perform poorly. Thus, we only use this approximation 
at distances where the error bound is small. All unit cells that are close to the point 
of evaluation will be calculated following the exact explicit sum. In the following 
sections we will calculate the electric field from all unit cells in 2. Additionally, we 
will derive a bound on the error in the electric field from all unit cells. We will use 
that error bound to determine where our approximated method can be used in {2 for 


a specific evaluation point «°. 


4.3.5 Contribution from All Unit Cells in a Body 


Working only with the dipole term to start, we want to calculate the total electric field 
from all unit cells located at a’ € 2. Due to the fact that the governing electrostatic 


equations are linear, we can simply add the contribution from each individual unit 
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cell to get the total electric field. 


Arey Berra?) x S- Arey Ee" (x?) 


w/EQ 
2 (2? = x!) (2? xt) — 776i; 
ae ! 
S- p(x )l |x? _ a! | (4.19) 
ZED 


We pass from the sum to an integral over 2 by assuming that || > 1°. Additionally, 


we convert the integral to its div p form using integration by parts. 


Q ae — 2x 
B ! B ! 
| ee | ee a Aang | ee: 
=f V 5pj(x') (335) d°ax +f mle )nj (3 =.) ds 


We now work with the higher order terms, which we call the error, W(x"), and can 


write as a series of terms. 


2 |a? — x’ | 


W (a?) = | AG) =GGelViV;Vi (az) +e. 
Q 
2 


1 
4 _ * —— el 
UViV5 Vi & =) +... (4.21) 


Q= | le(y)ld (4.22) 


By summing over the error from all unit cells in (2 and again passing from the 


summation to an integral by assuming || > 1°, we arrive at the following. 


Of AES 
U (a) < = | Sos (4.23) 
amy |e? — a /P 
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We note that in the limit of | — 0 all error terms go exactly to zero. The approxi- 
mation becomes exact in the limit as expected and shown previously in [35, 49]. As 
|a — x’|? gets larger, the relative error decreases. We capture this relationship by 
breaking up the integral over (2 into n shells of unit cells. In each integral, we will 
assume a fixed value of n for the shell and the maximum error as a function of nl. 
nl in this case, is the smallest scaled value of |a’ — a’|? in a given shell, which leads 


to the largest possible error. We will also assume for the time being that 2 is in the 


entire space of R°, thus we will sum over n shells from a minimum nmin to oo. In 


the case of a real system, the summation can simply be cutoff when the n™ shell is 


outside of the actual domain 2. 
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as o (((n + 1)l)* — (nl)*) (4.24) 


(((n + 1)l)* — (nl)’) 


dr) a ; 
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Now if we take the limit of n — oo, we obtain the same result of the error being 
identically zero as when the limit of 1 + 0. We have just slightly tweaked our 
perspective. Now instead of shrinking the unit cell size, we are increasing the distance 


|? — a’| between the unit cell location and the evaluation point. The electric field 
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can now be written in terms of the div p term along with an associated error. 


Egprree (98) ~~ | V; (a) ee 7 vi da! 
: an iP Aréo|x® — x! 
+ | (a!)n,; me \ ag (4.25) 
oe Pj j Arey|x8 — x’) : 
ae on + an 
B —= 
U(x") <C, 5° 5° (4.26) 


Nmin P=4 


This approximation will only be used in regions far enough away from the evaluation 
point to maintain the total error in the electric field below an acceptable tolerance 
level. 

We break up 2 into 3 distinct domains for each evaluation point, x, see Figure 
4.3. Qatomistic i8 the region where the QC meshes are fully resolved and defects or 
external loadings exist. For example, if a charge vacancy is modeled, the surrounding 
region would be fully resolved. We note that in the fully resolved region, we relax 
the assumption of charge neutrality because the full summation is used and the 
polarization does not need to be defined or calculated. The electric field from this 
portion of (2 will be included using the exact summation. S2,ne is a set of shells of 
unit cells that are centered at x°. The size of this region is determined by the input 
acceptable error tolerance. All unit cells within this domain will be included exactly 


as well. 


: Bgl ayy, 
Bett (gy) — | | p(x wie = us) By Aap (4.27) 
QatomisticU2shell Areo|ax = de = y| 


The contribution to the electric field from all other unit cells in (2¢g, will be approx- 
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Figure 4.3: Full (2 domain broken down into fully resolved region Qatomistic, & Shell 
region (Qsney centered at a point, x, where the field is to be evaluated, and the far 
field region Qy,. A single unit cell located at a’ in the far field region is also shown. 


imated using Equation 4.25. 


Qfar = 2 \ Qatomistic U Vsheti (4.28) 


The full electric field can now be written as the combination of the contributions 
from the approximated far field domain and the exact atomistic domains with the 


associated error in Equation 4.26. 
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We note that the only source of error is from the approximated region, which is 
controlled by the size of the shell region. The full atomistic regions contribute zero 
error to the electric field, because they are exactly summed. In section 4.5, we 
will evaluate the error both analytically and numerically in the electric field using 


material properties for GaN at a variety of different n shell region sizes. 


4.4 Numerical Implementation 


A detailed description of the QC code can be found in [40, 50], we will just highlight 
the basics of the code and the additions necessary to calculate the electrostatic forces. 
The non electrostatic portion of the QC code is a mix of C, C++, and Fortran, with 
a majority of the code written in C. The code takes lattices as an input and meshes 
each lattice separately. Data structures including hash tables and neighbor lists are 
used to store the required information for calculating forces on atoms. A variety of 
mechanical loads can be applied and the code solves the forces equal to 0 in the body 
subject to these external loads. 

We implemented our electrostatics method in the code by creating a C++ sin- 
gleton class that handles all electrostatic calculations. The base QC code then calls 
methods that update the electrostatic data whenever the electric field is required. 
The electric field is calculated in three steps upon request from the base QC code. 
The first step is to account for the contribution from all exact summations, see Equa- 
tion 4.27. This region includes shells of unit cells around an atom of interest and 
all atomistically resolved domains in the triangulated meshes. We use hash tables 
keyed by the triple integer of lattice coordinates unique to each atom to store lists of 
atoms in these region for quicker evaluation. The size of the shell region is controlled 
via an input integer that denotes the number of shells of unit cells to include. The 
contribution to the electric field from this exact region is then summed using the 
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standard kernel. 

The second step is to include all contributions to the electric field from the ap- 
proximated far field interactions. Polarization density is calculated using Equation 
4.17 throughout 2. We assume constant polarization in each element of a triangu- 
lated mesh. We then use the polarization density field to calculate the electric field 
with the first part of Equation 4.25. The far field contribution thus evolves into 
calculating the jump in polarization, [p2— pi] = 0, between elements and the diver- 
gence on all boundary terms, p-n. The exterior boundaries are simply controlled 
with input options of a free surface boundary, a charge-neutralized (infinite domain) 
boundary, or a specified charge density on the boundary. The contribution to the 
electric field from the jump in polarization terms and the exterior boundaries are 
solved using Equation 4.25 with a basic adaptive numerical integration scheme. The 
interior boundary between the approximated and the exact region, however, requires 
additional work. 

In Section 4.3.5, we took the full summation of shells of unit cells around the 
atom of interest, essentially in a spherical manner. However, implementation of this 
spherical summation of atoms is not truly spherical when interfacing with Cartesian 
based unit cells. The boundary of the full atomistic region approximates a sphere, 
but is not exact. As the number of shells included in the full atomistic region is 
increased, the boundary approximation becomes better. For the range of shells that 
we are interested in, however, this approximation is still not good enough and the 
additional numerical error is large. In principle, this issue shouldn’t exist, because 
the far field terms would approximate the correct boundary in conjunction with 
the exact region. Implementation proved too challenging to eliminate the excess 
numerical error, though. Instead we use trapezoidal shells of unit cells, which are 


much simpler to implement and do not have the additional numerical error. The 
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additional numerical error does not exist, because we precisely match the shape of 
the interior boundary in the far field region with the exact shell region. 

The far field approximation for the interior boundary amounts to finding the inter- 
section of a trapezoidal box with a triangulated mesh of potentially rapidly changing 
element size. A p-n contribution to the electric field then occurs on this intersection, 
where the polarization comes from the intersected element. We solve this intersec- 
tion for all required field calculations with the use of CGAL|1], a computational 
geometry library. The contribution to the electric field is then computed using the 
integral in Equation 4.25 with a basic adaptive numerical integration scheme. The 
CGAL library unfortunately introduces several additional library dependencies into 
the code, but gives the advantage of greatly speeding up the calculation of the trape- 
zoidal box versus triangulated mesh intersections. The full exact region geometry, 
does not match our original derivation, but we are including more unit cells in this 
region and less in the approximated region. Thus, the only consequence is that our 
error bound is higher than necessary. 

The final step is to calculate the contribution to the electric field from external 
loadings. These loading can come in the form of an applied field on the body, a 
set of point charges near a surface, or a vacancy/void in the body. In the case of 
an applied field the proper terms are directly added, while the standard kernel is 
used to add point charge interactions to atoms in 2. The void/vacancy case is not 
really an external loading, but is a defect which is natively simulated in the fully 
resolved portion of the mesh. Because this defect is in the exact summation region, 
no changes are required, the atom or atoms in question are simply not included in 
the summation. 

We note that the above steps are done for every representative atom in (2. If 


the electric field is needed at an atom, x°, that it is not a representative atom, 
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the electric field is interpolated at that location using shape functions. We assume 
that the interpolation of the electric field is valid and expect this assumption to 
be acceptable. We note that at any locations where the electric field is needed, 
the deformation field and hence the electric field will be slowly varying. With the 
electric field calculated , the base QC code requests the electric field at all required 
atoms and calculates the force by multiplying by the charge on the atom. This 
force is then added in the relevant places of the base QC code to the short-range 
forces from the other portion of a given interatomic potential, see Equation 4.8. The 
total resulting forces are then minimized to less than a given tolerance level with a 
conjugate-gradient algorithm. 

All electrostatic calculation in the code are implemented with threading for high 
performance computing. We use a queue structure for each thread to grab addi- 
tional chunks of work. There are multiple potential optimizations in terms of work 
chunk size, the number of threads, and the underlying data structures sizes. The 
optimal parameterization will depend on material properties, system geometry, and 


computational hardware, among other factors. 


4.5 Results 


We now present convergence results for our electrostatic summation method. All 
results are for a single crystal of wurtzite GaN. The wurtzite crystal structure has 
4 atoms per unit cell, 2 Ga and 2 N atoms. We assume a crystal structure parame- 
terization from a GaN core-shell potential [83]. GaN is a HI-IV type semiconductor 
with polarization in the 0001 direction [77]. The effective charges are +2 for Ga, 
-2.5 for the N (electron) shell, and 0.5 for the N core. All of the following results 
assume charge neutralized boundary conditions, thus the size of the simulation does 
not matter, only the size of the exact summation region. 
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We start by computing the formal error bound from Equation 4.26. The C 
constant is calculated from the GaN wurtzite structure. We calculate the total error 
in the electric field at a single evaluation point assuming an infinite surrounding 
crystal. The total error converges in Figure 4.4 to a singular value as expected from 


the convergence properties of Equation 4.26. We also calculate the error bound if 


40° Convergence of Total Error 


Total Error of Nth Shell and Lower 
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Figure 4.4: Convergence of total possible error in electric field from error bound in 
Equation 4.26 for an infinite crystal. Each successive shell adds less total error as 
expected because of the increased distance from evaluation point. 


a calculation is performed with the exact atomistic region the size of N shells of 
unit cells in Figure 4.5. Hence, if le~? is an acceptable error bound on the electric 
field, around 13 shells of unit cells surrounding the evaluation point need to be 
exactly summed. The contribution to the electric field from the rest of 2 can use 
the approximated method. 

We now calculate the error in the electric field from the actual implementation in 
the QC code. GaN in the wurtzite configuration only has polarization in the 0001- 
direction, so the interior boundary with the exact atomistic region will only have 
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Figure 4.5: Error bound for electric field with an exact atomistic summation per- 
formed with N unit cells. 


a charge distribution on the 0001 faces of the trapezoidal shell. Figure 4.6 shows 
the charge density from the interior boundary of the exact atomistic region and the 
approximated region centered about an evaluation point a. Figure 4.7 shows the 
contribution to the electric field in the 0001-direction from the same interior bound- 
ary. Figure 4.8 and Figure 4.9 are also provided to spatially show the contribution to 
the electric field from the far field region with different levels of integration. Figure 
4.10 shows the convergence of the electric field calculations for a single evaluation 
point. The three line plots show the different levels of adaptive quadrature used to 
solve Equation 4.25. The resulting convergence is well below the error bound from 
Figure 4.5 and decreases as expected. We also show the convergence of the exact 
shell region. It stands to reason that if the far field contributions converge to a value 
that the exact shell region must as well. This convergence is shown in Figure 4.11 


and is very similar as expected to the far field region convergence. 
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4.6 Discussion 


We have developed a method for calculating electrostatic interactions without as- 
suming periodic boundary conditions. Our method is a combination of FMM and 
the RF method. The RF method assumes a homogeneous domain outside of an ex- 
actly computed region. We do essentially the same thing, except we do not assume 
a homogeneous domain and instead allow the body to vary accordingly with the 
external loadings. The far field terms are calculated using only the dipole moment, 
which is essentially a simple form of FMM. The method has been implemented in a 
high performance computing framework and coupled with an existing QC implemen- 
tation. Currently, our electrostatic method has only been implemented with the one 
QC code. However, our electrostatic method is compatible with all QC methods and 
can be implemented in other methods outside of the QC specific field. A variety of 
electrostatic and mechanical boundary conditions can be applied to a system. Ad- 
ditionally, the error in the electrostatic method can be controlled by the size of the 
exact atomistic region surrounding a point of evaluation. We showed the resulting 
error both analytically and numerically. The numerical error is below the analytical 
error bound as required. 

This method for performing electrostatic calculations provides a high degree of 
control to a user. Users need to weigh the computational cost, numerical error, 
method error, problem size, error tolerance, and computational hardware when run- 
ning simulations. A number of scenarios that a user might encounter are highlighted 


below. 


e + exact atomistic region, | method error, + computational expense 
e ¢ numerical quadrature, | the numerical error, } computational expense 


e + defect size/effective range, + problem size, small + method error, + computa- 
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tional expense 


A number of other scenarios exist, but the end result is that a user has to make 
the decision for a variety of input parameters depending on the specific problem of 
interest. We are currently using the method to apply various types of loadings on 


an ionic solid and are investigating the computational scaling of the method. 
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Charge Density on Interior Boundary 


40 0.0 


Figure 4.6: Charge density on top 0001-plane of interior trapezoidal boundary be- 
tween the exact atomistic region and the approximated region. The size of the exact 
atomistic region is 8 shells of trapezoidal unit cells. 
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Contribution to Electric Field from Interior Boundary 
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Figure 4.7: Contribution to electric field in 0001-direction from top 0001-plane of 
interior trapezoidal boundary between the exact atomistic region and the approxi- 
mated region. The size of the exact atomistic region is 8 shells of trapezoidal unit 
cells and the integration level is 1. 
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Contribution to Electric Field on Interior Boundary 
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Figure 4.8: Contribution to electric field in 0001-direction from top 0001-plane of 
interior trapezoidal boundary between the exact atomistic region and the approxi- 
mated region. The size of the exact atomistic region is 8 shells of trapezoidal unit 
cells and the integration level is 2. 
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Contribution to Electric Field on Interior Boundary 
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Figure 4.9: Contribution to electric field in 0001-direction from top 0001-plane of 
interior trapezoidal boundary between the exact atomistic region and the approxi- 
mated region. The size of the exact atomistic region is 8 shells of trapezoidal unit 
cells and the integration level is 3. 
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Convergence of Electric Field 
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Figure 4.10: Numerical convergence of electric field contribution in 0001-direction 
from far field region using QC code with an infinite single crystal of wurtzite GaN. 
Integration notation denotes an increased used of numerical quadrature points used 
in solving the contribution from field terms. 
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Convergence of Electric Field 
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Figure 4.11: Numerical convergence of electric field contribution in 0001-direction 
from exact shell region using QC code with an infinite single crystal of wurtzite GaN. 
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Conclusions and Discussion 


Our focus in this thesis was to extend multiscale methods, specifically the QC method 
to long-range electrostatic interactions without assuming periodicity in ionic solids. 
We achieved this in Chapter 2 by using polarization as a multiscale linking quantity 
between atomistic and continuum scale electrical quantities. Our coarse-graining 
strategy relies on ideas developed by [35] and others that followed them [63, 80, 59]. 
The method that we have presented enables QC and other multiscale methods to 
go beyond purely short-range interactions, that are characteristic of many struc- 
tural materials, to functional and electronic materials of interest today. We used the 
method to investigate some electromechanical properties of an ionic solid. Addition- 
ally, we discuss the importance of boundary conditions and the details of calculating 
the polarization. 

In Chapter 3 we extend a non-local cluster force QC method originally developed 
by [40] to multiple lattices. We did this by treating the problem as the union of many 
simple lattice problems, some of the work follows unpublished results from [42]. We 


test and verify the implementation in a high performance computing code. Results 


118 


from the relaxation of nanovoids from a vacancy to 1nm are presented. These voids 
are also relaxed subject to various increments of volumetric expansion. We show the 
centrosymmetry parameter changing between the relaxed initial configuration and 
the relaxed configuration subjected to volumetrically expanded boundary conditions. 

In Chapter 4 we extend the results from Chapter 3 to ionic solids. We add 
electrostatic interactions by calculating the full charge-charge interaction summation 
for atoms that are close to the evaluation point. Farther away, we develop and use a 
method following ideas from the FMM [28, 7, 14]. Essentially, we include interactions 
from the approximated region with just the first term (dipole) from the multipole 
expansion. All higher order moments are treated as a bound on the maximum error 
in the approximated region. We then calculate the error bound both analytically 
and numerically for GaN and compare the results. 

In the Appendices, we write a manual for using the simulation code from Chap- 
ters 3 and 4. We are continuing to use the code to investigate interesting physical 
phenomenon in a variety of materials and hope that future users of the code find it 
useful. This thesis also motivates future work including, tighter error bound deriva- 
tions, computational scaling tests, and inclusion of higher order moments in the 


approximated region to reduce the size of the exact region. 
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Appendix A 


‘Transformation of Unit Cell in a Crystal Lattice to 
Obtain a Unit Cell with One Face Parallel to a 


Given Plane 


We show in this section that any crystal lattice can be described by a unit cell with 
one face parallel to a given plane. This is an essential ingredient of our proof that 
the change in polarization due a change in the lattice unit cell is balanced by a 
corresponding change in surface charge density. 


We first show the following proposition. 


Proposition A.0.1. Consider a crystal described by lattice vectors {f1, fo, f3} and 
reciprocal lattice vectors {f', f?, f?}. Consider a rational plane’ with normal n. 


This lattice can also be described by lattice vectors {gi,g2,g3}, where go 1 n, 


' A rational plane has a normal n that can be represented n = 5°, nf’ with n; integers. 
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g3 = fz, and g, will be shown below to exist. 


Proof. Result 3.1 of [9] states that a crystal described by lattice vectors {f1, fo, fs} 
can equivalently be described by lattice vectors g, go, g3 if and only if these satisfy 
=>, j ul f; with ul being a 3 x 3 matrix of integers with determinant +1. 

We have g3 = f3, therefore 3 = w3 = 0 and w3 = 1. Set go = Mi fi + Mofo with 
My = ae Cor and M2 = FEET This choice ensures the following: gz Ln, M, 
and My) are integers, and gcd(M,, Mz) = 1. It implies that wi = My, u3 = Mo, 3 = 0. 

The remaining step to complete the proof is to show that we can find g; which 
satisfies g; = >> j ul f; where yu! has the restrictions mentioned above. Write gi = 
Ni fi+ Nofot+ N3f3. Then det w = 1 gives us the condition NjMz— M,N2 = 1. The 
extended Euclidean algorithm provides the existence of integer solutions N; and Np» 
to this equation [44]. In general, this algorithm provides integer solutions m,n to 
the equation am + bn = gced(a, b) where a,b are given integers. While the algorithm 
is constructive and ensures existence of solutions, it does not provide explicit forms 


for the solutions. 


Therefore, there exists g; = Ni f; + Nof2+N3f3 where N,, No are obtained from 


the algorithm above and N3 can be arbitrary. 


We wish to show that given a crystal described by lattice vectors { fi, fo, f3}, we 


can equivalently describe the crystal by lattice vectors {hi, ho, hg} with hy L n and 
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h3 | n. The vector n is the unit normal to a rational plane. 

We simply apply the proposition above twice in succession. First, we use the 
proposition to transform from { fi, fo, f3} to {gi, g2 1 n,g3 = f3}. We then use the 
proposition again to transform from {g1, g2,g3} to {hi, he = go,h3 Ln}. Therefore, 
the final description has {h1,h). L n,h3 | n}. The condition det 4 = +1 ensures 


linear independence; in fact, it preserves the volume of the unit cell. 
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Appendix B 


How ‘To Run QC Code 


In this appendix, we will provide documentation for running the QC code. We 
start with the contents of the main folder, B.1. Inside this folder are the main 
QC source files, a directory of input files with inputs for NiAl and GaN, a build 
script for the source files, and a .tar.gz of everything for easy movement. All of 
the files for compilation are in the quasi-src folder, B.2._ The sre sub folder is a 
static library of QC functions with all header files in the include folder. The mesh 
generation folder does exactly as it says and has a variety of mesh generation files. 
The triangulation_server folder serves as a server for building meshes in the main 
code. The applications folder contains the main driver folders and files for the QC 
program. The file CGAL-4.4.tar.gz is the repository for CGAL, directly downloaded 
from the CGAL website [1]. All other folders contain information and files, but are 
not needed by a typical user. 

The first step of using the QC code is to build it. We have made a script for 
automatically taking care of the build process for several architectures, B.3 and B.4. 


In the first part of the script, the user needs to specify the absolute path of the 
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quasi-v1.1.1 


<€ | > feHome software Quasi quasi-v1.1.1 quasi-src 


Places — my 
© Recent ual alll : _— 
f Home input_files build-quasi.sh quasi-v1.1.1.tar.gz 
[mi Desktop 
DB Documents 
~ Downloads 
dd Music 
©) Pictures 
H videos 
@ Trash 

Devices 


(a) Floppy Disk 

() computer 
Network 

@? Browse Network 

D connect to Server 


|_“quasi-src” selected (containing 28 items) _ 


B.1: Main QC folder. 


main quasi folder with quasi_loc. The build directory also needs to be specified with 
an absolute or a relative path by defining build_directory. Additionally paths for 
dependencies and libraries need to be included with -I or -L options, if they are not 
available system wide. This is shown for several architectures including a standard 
desktop with Ubuntu and Blacklight, an XSEDE resource. If the code is being run 
on a machine with modules then the modules just need to be added. The code has 
the following dependencies mpfr, qt (>4.0), gmp, boost(>1.50), cmake (>3.0), and 
gcc (>4.8) Lastly, the Fortran, C, and C++ compilers that are to be used need to be 
specified. The second part of the script takes care of the build process for CGAL, the 
static quasi library, and the main driver code under applications/nano-indentation in 
the build directoy. Various levels of output can be shown in the terminal or output 
to files in the relevant build director. The goal of this script file is to ease the build 
difficulty, though it is not a black box. Users need to be somewhat fluent in linux 
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quasi-src 


fe Home software Quasi quasi-v1.1.1 quasi-src 
Places — — = 
© recent ould al aul | 
7 fHome applications autom4te.cache bin get_defects 
[mi Desktop = = - = 
| Dioocuments nlf al aul mall 
a ~ Downloads include input_files lib material_data 
dd Music - = 
arPcares wal wal wal wl 
H videos mesh_generation new_mesh_ solver src 
@ trash generation 
pee -" -" mal | 
(a) Floppy Disk 
a system triangulation_ utils aclocal.m4 
Computer server 
Network a FP 
—_o — 
@? Browse Network aS — = 1 
unde 
D connect to Server AUTHORS CGAL-4.4.tar.gz ChangeLog config.h.in 
dnt & Basic 
dni T sssss 
= Coov | Th 
configure configure.ac COPYING INSTALL 
## To # Mak Direc 
biel # @co === 
ACLOC 
SUBDI # Coo | | Backu 
Makefile.am Makefile.in NEWS README 


B.2: Source file folder. 


commands to build the program. 

With the code built, several input files need to be generated before simulations 
can occur. The first step is to generate the mesh files. We have solely been using 
the indent_fcc file in the main input file folders, B.5, and subsequent code in quasi- 
src/mesh-generation, though others could be used. This input and code needs to be 


run once for ever lattice in a given system. We show a typical one for GaN. 


e material: specifies the lattice type, in the GaN case it would be Ga, N core, 


or N shell (actual name does not matter, but helps for clarity). 


e atomistic: denotes a region that is fully resolved a number of shells away. In 
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it! /bin/bash 


# set options 

quasi_loc="insert src location" 
#quasi_loc="/usr/users/8/jmarsha2/quasiSRC-5-28-14" 

numThreads="1" 

runCommand="./quasi" 

build_directory="./quasi-v1.1-build" 

# dependencies boost, mpfr, qt4, gmp, cmake 
#directory_for_dependencies="-I/usr/local/lib/" 
directory_for_dependencies="-I/usr/local/packages/Boost/boost_1_50_0/include" 
directory_for_libraries="-L/usr/local/packages/Boost/boost_1_50_0/1ib" 
#PATH=/data/myscripts:$PATH 

#export PATH=$PATH: /usr/local/packages/Boost/boost_1_50_0/1lib 

# 

# gcc compilers defined in shell 

# 

#Fortran_Compiler="gfortran" 

#C_Compiler="gcc" 

#CXX_Compiler="g++" 


# 

# gcc compilers built locally (Desktop) 

an 

# dependencies boost, mpfr, qt4, gmp, cmake 

#Fortran_Compiler="gfortran" 

#C_Compiler="gcc" 

#CXX_Compiler="g++" 
#Fortran_Compiler="/home/jason/software/gcc4.9.OBUILD/bin/x86_64-unknown-linux-gnu-gfortran" 
#C_Compiler="/home/jason/software/gcc4.9.@BUILD/bin/x86_64-unknown-linux-gnu-gcc" 
#CXX_Compiler="/home/jason/software/gcc4.9.@BUILD/bin/x86_64-unknown- linux-gnu-g++" 


# 

# gcc compilers built locally (Blacklight) 

# 
#Fortran_Compiler="/usr/users/8/jmarsha2/gcc4.9.0/bin/x86_64-unknown-linux-gnu-gfortran"” 
#C_Compiler="/usr/users/8/jmarsha2/gcc4.9.0/bin/x86_64-unknown- lLinux-gnu-gcc" 
#CXX_Compiler="/usr/users/8/jmarsha2/gcc4.9.0/bin/x86_64-unknown- Linux-gnu-g++" 

#export LD_LIBRARY_PATH=/usr/users/8/jmarsha2/gcc4.9.0/1ib64 


# 

# gcc compilers in module (Blacklight) 
a 

source /usr/share/modules/init/bash 
module load mpfr/3.1.2 

module load qt/4.6.3 

module load gmp/6.0.0 

module load boost/boost_1.50.0 
module load cmake/3.0.1 

module load gcc/4.8.2 
Fortran_Compiler="gfortran" 
C_Compiler="gcc" 

CXX_Compiler="g++" 


B.3: Build script, part 1. 


our example, we specify a 2x2x2 unit cell radius, so the total atomistic region 
would be 5x5x5, one unit cell as the base, and two in the plus and minus 


directions. 
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HHHHHHHHHHAHAHHHAAAA AAA RRRR RRR 
# SHOULD NOT NEED TO MODIFY BELOW HERE 
HHHHHHHHHHHHHHAHAAAAAAAAA ARRAS 


# get pwd and date 
MY_PWD=$(pwd) 
dateDay= date +%Y-%m-%d~ 


# build CGAL 


echo VHKEKKKKKEKEKKEKEEEEKEKEEEKEEEREEEKEKEKEKEKEKEREEKKEE ' 


echo ' Building CGAL 4.4 ... ; 
echo VAKEKKKEKKEKEEEKEREKKEKEEEEEEKEKEEEKEEEKEKREEEKEREEKEEE ' 
cd Squasi_loc 

rm -rf CGAL-4.4 

tar -zxf "CGAL-4.4.tar.gz" 

cd CGAL-4.4 
cgalconfigOUT="cgal-configure-results-$dateDay. txt” 
cgalmakeOUT="cgal-make-results-$dateDay. txt" 

cmake -DBUILD_SHARED_LIBS=FALSE . #> $cgalconfigOUT 2>&1 
make #> $cgalmakeOUT 2>&1 


echo VHKKKKKEKKEKKEEEEREKKEKEKKEEEEKEEEKEEEEEEKEKEEKEKREEKKEKE ' 


echo ' Building CGAL 4.4 ... Done 1 


echo VHKEKEKEEKKKEEKEKEKEEKEEKEEKEEEEKKEEKEKEEEEEKEEKEKEREKEEKEE | 


# build directory 

cd SMY_PWD 

mkdir $build_directory 
cd S$build_directory 


quasi_configure=$quasi_loc"/configure" 

C_FLAGS="" 

CXX_FLAGS="" 

CPP_FLAGS="-std=c++11 -02 -frounding-math "S$directory_for_dependencies 
LD_FLAGS=S$directory_for_libraries 
configOUT="quasi-configure-results-$dateDay. txt" 
makeOUT="quasi-lLib-make-results-S$dateDay.txt" 
makeOUT2="quasi-main-make-results-S$dateDay.txt" 
outputOUT="results-$dateDay.txt" 


# build quasi 


echo VHEKEKKKEKEKEREKEREREREREEREREREREREREREREREREREEKE | 


echo ' Building Quasi 1.0... ‘ 

echo VHKEKEEEKEKEKEEEEEKEEEEEEEEEKEEEEEEEEEEEKEKEEEEKKEE ' 

F77=$Fortran_Compiler CC=$C_Compiler CXX=$CXX_Compiler CFLAGS=$C_FLAGS CXXFLAGS=$CXX_FLAGS CPPFLAGS= 
SCPP_FLAGS LDFLAGS=$LD_FLAGS $quasi_configure #> $configOUT 2>&1 

make -j SnumThreads #> $makeOUT 2>&1 

cd applications/nano-indentation 

make -j SnumThreads #> $makeOUT2 2>&1 


echo VHKKEKEEEEKEEEKEEEEKEEEKEEEEEEEKEEEEEKEREREEEEEEEEKE | 


echo ' Building Quasi 1.0 ... Done : 
echo VHEKEKEKEKEKEKEREREREREREEREEEEKEEKREREREREKEREREREEKE ' 


B.4: Build script, part 2. 


e full: specifies the full system size and will be calculated below. 


e indentor: specifies the location, if the indentor is to be used, if not any values 


are acceptable. 


e lattice_type: exactly what you would expect. 0=fcc, 1=bcec, 2=wurtzite, oth- 


ers can be added in the code. 
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# GaN 


material Nc 

atomistic 2.822 

# zis half, if problem_size == 1 
full 24 24 12 
indentor 00 

output GaN_6.inp.gz 
Lattice_type 2 


Lattice_constant 3.229234 3.229234 5.162889 


#sum of n_levels and atomistic * 2 must equal full atomistic 
#and must be divisible by last level size 


mesh_levels 2 
n_level 0 4 
n_level i: © 
#n_level 24 
#n_level 3 16 
#n_level 4 48 
#n_level 5 48 
#n_level 6 48 
#n_level 7 48 
#n_level 8 48 
# problem size 0 = non-mirrored z plane, 1 = mirrored z plane 
problem_size 1 


# zplane © = normal creation, 1 = only add nodes in z plane 
zplane 0 


B.5: Mesh generation input file. 


e lattice_constant: relevant lattice parameters. 
e mesh_levels: is the number of coarse graining that are to occur. 


e n_level: two parameters must be specified, the level number and the jump in 
unit cells. The full size must be calculated from these levels and is twice the 
atomistic plus the sum of the level jumps. Hence 2 x (2+ 4+ 6) = 24, so the 
full size would be 24x24x24. 


e problem size: specifies if the full mesh is to be used or only a half mesh. 


Hence, if only a half size is used (0), the full problem size is now 24x24x12. 
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Note that the full domain size must be divisible by the last coarse-grain jump 


size. 


e zplane: is used to build pillars. The domain is essentially the size of the x and 
y atomistic region z tall. The full parameter would thus be the atomistic size 
in the x and y direction and the z direction comes from the number of specified 


coarse-grain levels. 


Additionally, the output parameter can be anything, but must end in .inp.gz. There 
needs to be one mesh file per lattice with the naming convention of ” name” .inp.gz for 
the base lattice and ”name” _n.inp.gz, where n starts at 2. These mesh files need to be 
moved to the main driver folder in the build directory, in this case applications/nano- 
indentation. 


We also need to make the input for the main driver code, B.6, B.7, and B.8. 
e number_threads: number of pthreads to solve the problem with. 
e periodic: not used. 


e data file: mesh input file, only list the base lattice input file, code will auto- 


matically read in all others assuming proper naming. 
e materials_file: not used. 


e indentor parameters: self explanatory. Specifics can be found in the indent.c 


source file. 


e void parameters: self explanatory. Specific void types can be found in inter- 


nal_void.c. 
e numLoads: number of load increments to apply to system. 
e boundaryFlag: switches between applying the load on all atoms, just the 


boundaries, or the incremental deformation is turned off. 
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cs 

# setup input 

ca 

number_threads 1 

periodic 

data_file GaN.inp.gz 
materials_file materials.dat 


- 

# indenter info (0=off,1=on) 

z 

radius 700 

constant 100 

displacement 0.0 0.0 -0.1 

position 1000000.0 1000000.0 1000000.0 
indentEnable (0) 

#position s/2 s/2 s/1 


q 
# void info (0=false,1i=true) 

# 

voidEnable (0) 

#voidCenter 77.54 44.64 87.2 
#voidCenter 697.536 401.76 747.97 
voidCenter 348.5 200.9 372.24 
#voidcenter 34.3044 34.3044 34.3044 
voidNumParams 2 

voidParams 8.53 16.0 

voidType V 

#voidNumParams 
#voidParams 
#voidType 


ANF 


oa 
# loading info (flags: -i=off, 9=all atoms, 1=boundaries) 
= 


numLoads | 

boundaryFlag -1 

Fi 0.0000 0.0000 0.0000 

F2 0.0000 0.0000 0.0000 

F3 0.0000 0.0000 0.0000 

# charge number, charge increment, charge location 
#charge 0 -1.0 52.8 38.64 430.0 


# relax atomistic electro, used for surfaces (0 do not relax, 1 do relax) 
relax_electro_atomistic 0 


B.6: Main driver input file, part 1. 


e F1-F3: incremental deformation matrix to be applied to the system. 


e charge: incremental charge to be applied to the system. (charge num, strength, 
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numLattice 6 

shift 00.0 0.0 0.0 -1 

shift 1 1.614617 0.9321995596 2.5814445 -1 
shift 2 0.0 0.0 1.990655112 -1 

shift 3 1.614617 0.9321995596 4.572099612 -1 
shift 4 0.0 0.0 1.987572867 -1 

shift 5 1.614617 0.9321995596 4.569017367 -1 


potential GaNCoreShellCorrect 10.0 0 
max_num_shells 32 


# electrostaticBcs 

# © - charge neutralized, 1 - free surface, 2 - specify charge density 
# value only used for case 2 

xBeginElectroBcs 
xEndElectroBcs 
yBeginElectroBcs 
yEndElectroBcs 
zBeginElectroBcs 


10] 
0 
0 
0 
8 
zEndElectroBcS 0 


# remove residual forces 

# 0 - do not remove, 1 - remove 
# all non atomic atoms 
allResidual Z 

xBeginResidual 1 

xEndResidual 1 

yBeginResidual 1 

yEndResidual 1 

zBeginResidual 1 

zEndResidual 1 


# 

# CG info 

# 

tolerance 1.0e-2 
maxIterations ic} 
debugLevel 2 
LineTolerance 1.0e-6 


LineIterations 0 
remeshTolerance 10.0 
#remeshTolerance 0.002 


# 

# neighbor list flag if ® do not save locations, if 1 do save 
# 

crossNeighListFlag 1 


# set eam method (0 = from neighbor lists, 1 = from cluster calculation) 
eam_method (0) 


# 

# process restart file for defects (0 off, 1 on) 
# 

getDefectsFlag 0 

defectCenter 45.7392 45.7392 45.7392 
defectBox 5.0 5.0 5.0 

defectLat 2.8587 


B.7: Main driver input file, part 2. 
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# 


# add node options (flag 9 = use cutoff number, flag = 1 use remesh number) 


a 


addNodeF lag 1 
addNodeRestartNum 0 
addNodeCutoff 50 


B.8: Main driver input file, part 3. 


X, y, 2). 
relax_electro_atomistic: 0 do not relax surface and apply negative dead loads 


to correct for boundary interactions, 1 do relax surface forces. 
numLattice: number of lattices in the unit cell or system. 


shift: shift for each lattice in unit cell, last parameter determines is used for 


core-shell models. 


potential: specifies the potential to use and cutoff radius. See ComplexIn- 
put.cc, input.c, and PairPotentials.cc for more details and for implementation 


of new potentials. 


max_num_shells: used to index shells of atoms, the larger the number the 


more required memory. Program will exit if not enough shells allocated. 
electrostatic bcs: self explanatory 


residual forces: used to remove forces, if the input structure is not at equi- 


librium. 
conjugate gradient: see CGNonLinearSolver.cc. 


crossNeighListFlag: switches between different methods for computing neigh- 


bor lists. Suggest only using 1 at this time. 


eam_metho: switches between different methods for calculating eam densities, 


suggest only using 0 at this time. 
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e get defect parameters: used for post-processing only, outputs .3D, and .pdb 
full atomistic files from specified region with centrosymmetry parameter. Used 


to see dislocation nucleation and surfaces. 


e remesh parameters: used to specify the largest number of nodes added before 


moving on to next step or number of times remeshing can occur. 


This file must be named quasi.ini. The code can be run in a standard way ” ./quasi” 
from the build directory (application /nano-indentation). The results of the code are 
output into a numbered folder with 2 main file types: the first is a node*.plt.gz file 
that has deformations and other variable information on the mesh. This file is easily 
viewed in Visit [15]. The second file type is a restart*.inp.gz, which outputs a con- 
figuration that can be restarted from. Additionally, a couple of run time options are 
available for the code, the most used are -n number_threads and -r restart_file_ name. 

Lastly, support for this code can be found from Jason Marshall at his email 


address: jason.p.marshall at gmail.com. 
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Appendix C 


QC Code Computational Scaling 


Scaling results for the QC code are presented in this Appendix. All tests were 
performed on Blacklight a high performance computing resource at the Pittsburgh 
Supercomputing Center under XSEDE. The specs for Blacklight are quoted directly 
from the website. ” Blacklight is an SGI UV 1000cc-NUMA shared-memory system 
comprising 256 blades. Each blade holds 2 Intel Xeon X7560 (Nehalem) eight-core 
processors, for a total of 4096 cores across the whole machine. Each core has a clock 
rate of 2.27 GHz, supports two hardware threads and can perform 9 Gflops. The 
sixteen cores on each blade share 128 Gbytes of local memory. Users can run shared 
memory jobs of up to 16 Tbytes.” [13] 

We test the scaling results only for the electrostatic algorithms in the code. The 
first result is for just the exact calculations. These calculations include the time 
spent to calculate all atoms within a shell region and any atoms in the fully resolved 
atomistic region, see Figure C.1. The next set of results are just for the far field 
contributions to the electric field. Each of the three figures show results for different 


levels of numerical integration. As expected, more integration points increase the 
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computational time, see Figures C.2, C.3, and C.4. 

The results upon initial inspection are not promising. Upon closer inspection, 
however, there is reason to be optimistic. For small clusters, there is essentially 
no speed up with any number of threads. This is expected because of the minimal 
amount of work. As the number of clusters and amount of work is increased, some 
speedup is noticed with more threads. The uneven speedup can be most likely be 
attributed to the lack of optimization undertaken in the testing. Thread affinity 
was not set and is probably playing a large role in the poor scaling results. While 
Blacklight can be viewed essentially as a single computer, where all memory can 
be accessed from all processing units, the time taken to reach a memory location 
is not equal. Even on a single blade, the communication cost between a core on 
one Nehalem processor pays a time penalty of 1.3 to reach memory on the other 
processor as compared to no penalty to reach memory on its own processor. Thus, 
a time penalty is paid if the memory is not localized. The test for scaling did not 
take advantage of memory locality. Additionally, without setting the thread affinity, 
the first 2 threads were most likely placed on separate Nehalem cores, which could 
explain the increase in time when going from 1 thread to 2. Not setting thread 
affinity also allows threads to move around between different processing units, this 
migration thrashes the cache and also contributes to the memory access issues. 

The code implemented has the capability to set thread affinity already built in. 
Additionally, data locality in the code can be implemented to minimize the memory 
access times. The ideal optimization of the above factors will be both problem and 
hardware dependent. Testing is ongoing to determine the optimal configuration for 


Blacklight and in the future other computational resources. 
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C.1: Scaling results for exact cluster summation. 
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Scaling of Field Calculation, Integration = 1 
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C.2: Scaling results for far field electric field contribution with a numerical integration 
level of 1. 
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Scaling of Field Calculation, Integration = 2 
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C.3: Scaling results for far field electric field contribution with a numerical integration 
level of 2. 
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Scaling of Field Calculation, Integration = 3 
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C.4: Scaling results for far field electric field contribution with a numerical integration 
level of 3. 


ELECTROELASTICITY OF POLYMER NETWORKS 
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Asstract. A multiscale analysis of the electromechanical coupling in elastic dielectric is con- 
ducted, starting from the discrete monomer level through the polymer chain and up to the 
macroscopic level. Three models for the local relations between the molecular dipoles and the 
electric field that can fit a variety of dipolar monomers are considered. The entropy of the network 
is accounted for within the framework of statistical mechanics with appropriate kinematic and 
energetic constraints. At the macroscopic level closed form explicit expressions for the behaviors 
of amorphous dielectrics and isotropic polymer networks are determined. None of which admits 
the commonly assumed linear relation between the polarization and the electric field. The anal- 
ysis reveals the dependence of the macroscopic coupled behavior on three primary microscopic 
parameters: the model assumed for the local behavior, the intensity of the local dipole, and the 
length of the chain. We show how these parameters influence the directional distributions of 
the monomers and the hence the resulting overall response of the network. In particular, the 
dependences of the polarization and the polarization induced stress on the deformation of the 
dielectric are illustrated. More surprisingly, we also reveal a dependence of the stress on the 


electric field which stems from the kinematic constraint imposed on the chains. 


1. INTRODUCTION 


Electro-active polymers (EAPs) are materials that deform in response to electrostatic excita- 
tion. Due to their light weight, flexibility and availability these materials can be used in a wide 
variety of applications such as artificial muscles (Bar-Cohen, 2001), energy-harvesting devices 
(McKay et al., 2010), micropumps (Rudykh et al., 2012), and tunable wave guides (Shmuel 
et al., 2012), among others. The principle of the actuation is based on the attraction between 
two oppositely charged electrodes attached to the faces of a thin soft elastomer sheet. Due to 
Poisson’s effect, the sheet expands in the transverse direction. Toupin (1956), in his pioneering 
analysis, found that this electromechanical coupling is characterized by a quadratic dependence 


on the applied electric field, and this was later verified experimentally (Pelrine et al., 2000b; 
1 


Noy Cohen, Kaushik Dayal, and Gal deBotton 2 


Kofod et al., 2003; Wissler and Mazza, 2007). However, EAPs have a low energy-density in 
comparison with other actuators such as piezoelectrics and shape memory alloys (Bar-Cohen, 
2001) and their feasibility is limited due to the high electric fields (~ 100 MV/m) required for 
a meaningful actuation. The latter is a result of the relatively low ratio between the dielectric 
and the elastic moduli (Pelrine et al., 2000a). Specifically, common flexible polymers have low 
dielectric moduli while polymers with high dielectric moduli are usually stiff. Nevertheless, a 
few recent works suggest that this ratio may be improved. Huang et al. (2004) demonstrated ex- 
perimentally that organic composite EAPs experience more than 8% actuation strain in response 
to an activation field of 20 MV/m. The experimental work of Stoyanov et al. (2010) showed that 
the actuation can be dramatically improved by embedding conducting particles in a soft polymer. 
In parallel, theoretical works dealing with the enhancement of coupling in composites also hint 
at the possibility of improved actuation with an appropriate adjustment of their microstructure 
(e.g., deBotton et al., 2007; Xie et al., 2008; Rudykh et al., 2013; Tian et al., 2012; Galipeau 
and Ponte Castaneda, 2012; Lopez-Pamies, 2014; Spinelli et al., 2015, among others). These 
findings motivate an in-depth multiscale analysis of the electromechanical coupling in EAPs 
which is inherent from their microstructure. 

We note that the response of polymers to purely mechanical loadings was extensively inves- 
tigated at all length scales. A macroscopic level analysis and models describing the behavior 
of soft materials undergoing large deformation, such as polymers, were developed by Ogden 
(1997). The microscopic level analysis of Kuhn and Griin (1942) yielded a Langevin based 
constitutive relation and paved the way to various multiscale models such as the 3-chain model 
(Wang and Guth, 1952), the tetrahedral model (Flory and Rehner, 1943; Treloar, 1946) and 
the 8-chain model (Arruda and Boyce, 1993). The use of statistical mechanics in the analysis 
of mechanical systems was discussed by Su and Purohit (2012), and by Warner and Terentjev 
(2003) for polymer-based liquid crystal elastomer active materials. 

The electric response of dielectrics to electrostatic excitation was examined by Tiersten (1990) 
and Hutter et al. (2006), among others, macroscopically as well as through their microstructure. 
Starting with the examination of a single charge under an electric field, the relations between 
the different macroscopic electric quantities, such as the electric displacement, the electric field 
and the polarization, and the microscopic quantities, such as the free and bound charge densities 


and the dipoles were defined and analyzed. 
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A macroscopic level analysis of the coupled electromechanical response of dielectrics began 
with the work of Toupin (1956). Later on, an invariant-based representation of the constitutive 
behavior of EAPs was introduced by Dorfmann and Ogden (2005) and extended to anisotropic 
materials by Bustamante (2009). The possible influence of the deformation and its rate on the 
electromechanical coupling was investigated by Zhao and Suo (2008); Ask et al. (2012b,a); 
Jiménez and McMeeking (2013); Gei et al. (2013) and Deng et al. (2014). Initial multiscale 
physically motivated analyses of the electromechanical response were performed by Cohen and 
deBotton (2014), Cohen and deBotton (2015) and Cohen et al. (2015). 

In this work, a multiscale entropy-based analysis of the response of a dielectric to an elec- 
tromechanical loading is introduced. Employing the dipole models proposed by Stockmayer 
(1967) and Cohen and deBotton (2015), we start by analyzing both the mechanical and the 
electrical behavior of the dipolar monomers. Next, in order to determine the configuration 
of a polymer chain, we maximize its entropy and find the most probable distribution of the 
dipolar monomers composing it under given chain constraints. An integration over all dipolar 
monomers distributed according to the obtained conformation yields the entropy of the polymer 
chain. Lastly, assuming that the response of the polymer and its chains is strictly elastic, we em- 
ploy the first and the second laws of thermodynamics and integrate over all chains to determine 
the macroscopic response. 

We begin this work with a short theoretical background, followed by a thermodynamical 
analysis of the response of polymers in accordance with McMeeking et al. (2007); Tian et al. 
(2012); Cohen and deBotton (2014) and Cohen and deBotton (2015). In section 3 we analyze the 
microstructure of a polymer chain using statistical considerations through the entropy. Section 
4 deals with the behavior of amorphous dielectrics subjected to a purely electrical loading and 
a comparison with a well-known macroscopic model is made. Next, we consider polymer 
networks under electro-mechanical loadings, and determine the evolution of the microstructure 


and the macroscopic response. A few conclusions are gathered in section 6. 


2. THEORETICAL BACKGROUND 


Consider the deformation of a dielectric body from a referential configuration to a current one 
due to an electro-mechanical loading. The body occupies a region Vo C R? with a boundary 0Vo 
before the deformation and a region V C R? witha boundary OV at the current configuration. xo 


and x denote the locations of the material points before and after the deformation, respectively. 


Noy Cohen, Kaushik Dayal, and Gal deBotton 4 


The mapping of the material points from the reference to the current configuration is x = (xo) 


and the corresponding deformation gradient is 


F = Vx, (Xo), (1) 


where Vx, is the gradient operation with respect to the referential coordinate system. The ratio 
between the volumes of an infinitesimal element in the current and the reference configurations, 
J = det (F), is strictly positive. The velocity of the material points is v(x) and the spatial 
velocity gradient is 


L=V,v =FF|, (2) 


where the gradient Vx is carried out with respect to the current coordinate system. 
The body is subjected to an electric field E(x), satisfying Faraday’s law Vy XE = 0 throughout 
the entire space. Consequently, a scalar electric potential ¢ can be defined such that E = —V x@. 


The electric displacement field is 
D(x) = eE(x) + P(x), (3) 


where €o is the permittivity of vacuum and P(x) is the polarization, or the electric dipole-density. 
If a linear relation is assumed between the polarization and the electric field, then P = y Ee.) E 
where y is the susceptibility. In this case D = e9¢-E where €, = y +1 is the relative permittivity. 
We recall that in vacuum P = 0. In ideal dielectrics or in a continuum with no free charges the 


electric displacement field is governed by the equation 
Vx D=0. (4) 


The electrical boundary conditions are given in terms of either the electric potential or the 
charge per unit area p,, such that D - fi = —p,, where fi is the outward pointing unit normal 
to the boundary in the current configuration. Practically, in common EAPs settings, pq is the 
charge on the electrodes (Carpi et al., 2010). The mechanical boundary conditions are given in 
terms of the displacement or the mechanical traction t. Due to the presence of the electric field 
in the surrounding space, on the boundary the true or the current stress in the body o satisfies 


the condition [ao — a0” | = t, where 


o” = (E@B- 5 (B-E)1), (5) 
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is the Maxwell stress in vacuum and t is the mechanical traction. In the body the stress satisfies 
the equilibrium equation 


Vx :o =0. (6) 
Following McMeeking et al. (2007), the first law of thermodynamics states that 


dw d 


5 fue. mae a] 5 =n: BO Sa ee a [samaw, (7) 


Vo Yo 

where H and S are the electrical enthalpy-density and the entropy-density functions per unit refer- 
ential volume, respectively, W is the work of the external electrical and mechanical sources and T 
is the absolute temperature. We point out that the rate of the electrical enthalpy-density function 
is proportional to the rate of the internal energy-density function (McMeeking et al., 2007; Co- 
hen and deBotton, 2014) and stems from the application of the electric field, hence it is assumed 
that the variations of the total electrical enthalpy are negligible, i.e. df y) 7 (F, E = 0) dVo = 
(Treloar, 1975). The entropy-density function depends on the internal structure of the dielectric 
and evolves with the external loading. 


The rate of work of the external sources is 


dw _ tevda~ f po Pa (8) 


where we neglect body forces and assume no free charges in the material (Dorfmann and Ogden, 


2014). Eq. (8) can be written as (Tian et al., 2012) 


ral o” -E@P):LdV+ i (o —o0"):LaV- i pee iba EdV. 


R3\V 
(9) 
The rate of the electrical enthalpy is given via (Cohen and deBotton, 2014) 
d 1 (OH (FE) .7 OH(F,E) . 
— | H(B,E)d\Y = | - | ———— F : L + ——_ -E]dV, 10 
dt / aeen / J ( OF JE (10) 
Vo Vv 
and the rate of the entropy is 
d 1 (OS (F,E) _7 OS(F,E) . 
— F,E = {| ~— |———F : L+ ———_ -E 11 
5 [ S® ak i oF * OE ve “ 


Vv 
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Substituting Eqs. (9), (10) and (11) into Eq. (7) and collecting the coefficients of L and E yields 


_ oe eg OE) | ee aa 
o=[[e ao’ -E@P 5 ( AF T AF Jr):vav+ f wo oa’): Ldv 
R3/V 
il OS(F,E) OH (F,E) ; 


We follow Coleman and Noll (1963), and postulate that the latter holds for every admissible 


process. Therefore, 


o=a"+E@P+0o", (13) 

—_ 1 (OH(E,E) _OS(B,E 
one} (ED SEB) 0 

a 1 (.OS(F,E) 0H (EE) 
p=; (roe - Se. = 


In the surrounding space 0 = a”. 

The stress in Eq. (13) can be viewed as the sum of three contributions. For conciseness 
we will refer to the stress component a” as the mechanical stress and to the component 
E ® P as the polarization stress. The first term accounts for the evolution of the internal 
structure of the dielectric following the electro-mechanical loading. The second term, the 
polarization stress, stems from the application of an electric field over a polarizable medium. 
A similar decomposition was encountered in the phenomenological invariant-based framework 
of Dorfmann and Ogden (2005) and in the multi-scale analysis of dielectrics of Cohen and 
deBotton (2014), among others. The decomposition of the stress in these works was to a sum of 
(1) amechanical stress attributed to purely mechanical sources, (2) a polarization stress resulting 
from the applied electric field, and (3) the Maxwell stress in vacuum, where the latter two are 
similar to the second and the third terms in Eq. (13). However, it is important to note that any 
decomposition is essentially artificial, since from a practical point of view only the total stress 
can be measured in experiments. Eq. (15) reveals that the polarization depends on the variations 


of the microstructure caused by the external loadings. 


3. ENTROPY-DRIVEN ELECTROELASTICITY OF A POLYMER NETWORK 


Polymers are commonly used dielectric materials with a multiscale hierarchical structure 


of a large number of polymer chains, where each chain is a long string of repeating dipolar 
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monomers. The monomers can move or rotate relative to their neighbors thus providing the 
chains a freedom to deform (Flory, 1953). In this work the basic building block is the dipolar 
monomer which is treated as a small rigid rod and as a single dipole from a mechanical and an 


electrostatic viewpoint, respectively. 


3.1. A chain of monomers. We begin our analysis with a single chain with n dipolar monomers. 
The length of the line segment between the contact points of a monomer with its neighbors is /. 
The chain is subjected to an electric field E = EE and a mechanical deformation such that at the 
current configuration the vector connecting its two ends is r. We prescribe a coordinate system 


{f, Y, Z\, where Y and Z span the plane perpendicular to E. In this system a unit vector is 
& = cos OE + sind (cos oY +sin $2) , (16) 


where 0 < 6 < zis the angle between é and the electric field and 0 < ¢ < 27 is the angle of the 
projection on the plane perpendicular to E with Y. We also define the differential solid angle 
dl = sin@ dé d¢. 

The number of possible configurations of the chain is 


n! 


Q = ——_,, 
I]; (n!) 


(17) 


where n is the number of dipolar monomers aligned along é in the range 0) <9 <0%+4d0 
and 6” < ¢ < ¢” +dé. For convenience we define é m as the unit vector corresponding to 0 


and 6“. Consequently, the entropy of the chain is 


Sc=kin(Q)=k (: In(n) —n- »> n® In(n®) + »» n?) (18) 


i i 
where Stirling’s approximation is employed and k is the Boltzmann constant. The chain is 


subjected to three constraints: 


ye =n, (19) 


» in®é® =r, (20) 


and 
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5 aoe = He, (21) 


where h“ is the electrical enthalpy of a monomer directed along é = and Hc is the enthalpy of 
the chain. 

We assume that the polymer chain occupies the most probable configuration { n,n, a 
under the given constraints, and therefore we are interested in maximizing the entropy, 


Sc=k (i (Q) +a (>: n) — | +T> (>: nde” = | + ap nOpO — He}) (22) 


1 


1 


where a, tT and y are Lagrange multipliers (e.g., Davidson, 1962). The derivative of Sc with 


respect to n is 


— a (-In (n) +at7- 24 yh) if (23) 
from which we determine 
n = exp ( +7: Eg? + yh) ; (24) 
Upon substitution of the latter into Eq. (22), the maximum entropy that can be achieved by the 
system is 
Sc =k (nIn(n)-an-t-5-yHe). (25) 


We follow the works of Kuhn and Griin (1942) and Treloar (1975), among others, and assume 
that the polymer chains do not interact with one another. Consequently, in a volume element 
dVo, the total entropy-density and the total electrical enthalpy-density functions are obtained via 
s= a ya) “ and H = VG A i respectively. Differentiating Eq. (7) for the first law of 


thermodynamics with respect to the enthalpy of the k-th chain we obtain 


OH os 
aH® ~~ aH®” = 
C C 
from which we derive the relation 
1 
ee 27 
- ET’ (27) 


where Eq. (25) is used. 
From the constraint Eq. (19), (23) and (27) 


~ h 
(i) _ . — a — 
) n\’ = exp(a@) [oe (= é i dr =n, (28) 


1 
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where Eq. (27) is used and the summation is replaced by an integral over the orientations of the 
monomers. Therefore, 
n 
= —, 29 
exp (@) a (29) 
where 


h 


z= | ex(r-2- 7] dI, (30) 


is the partition function. Subsequently, we define the probability density function (PDF) that a 
monomer is in the direction € and has an electrical-enthalpy h as 


p (Bn) = e(r-2- 4). (31) 


and we note that { pd = 1, and Eq. (24) can be rewritten as n = np. An implicit equation 


from which the Lagrange multiplier t is computed follows from constraint Eq. (20), 


| epar= us (32) 
nl 


where again the summation is replaced by an integration over all possible orientations of the 
monomers. We note that n/ is the contour length of the chain. 

In order to compute the polarization and the total stress acting on a polymer according to 
Eqs. (15) and (13), respectively, we first formally differentiate the entropy of the chain Eq. (25) 


with respect to a given quantity e, 


Benn (-(2 OT “| t Or 1 am). (33) 


de je Oe i) f ae kT oe 


By differentiating Eq. (28) with respect to e, it follows that 


0 Oa OT ns 1 Oh 
nge(fvar)=n (52 fare. f evar oo 7a Pa] (34) 


oy ee r n Oh 
Oe Oe | kT Ha 


where Eq. (32) is used, and consequently 


n Oh Oa Ot Yr 
Vr) ae oe oe F G5) 
Substituting Eq. (35) into Eq. (33) yields 
OSc 1 (0He Oh kTt Or 
= -n | —pdr- cay 
de T | oe ee 7 | (6) 
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3.2. The polymer network. Following the premise that the chains do not interact with one 
another we have that a ag! 

de dV », je” @) 
where the summation is carried over all the chains. 


Making use of Eqs. (36) and (37), Eq. (14) can be written as 


(kK) (k) (k) 
kT 
ghee D(H sp par] iy A cs jer (38) 
k 


~ TdVo i OF 


The mechanical stress is composed of two parts: the first is related to the change in the 
electrical energy of the single monomer with respect to the mechanical deformation, and the 
second accounts for the mechanical forces deforming the end-to-end vector. Note that in many 
materials and regimes of practical importance, it is reasonable to assume that the monomer 
is rigid compared to the polymer chain. In that case, the enthalpy of the monomer does not 
depend on the deformation gradient, and the variation with respect to it vanishes. We also recall 
that previous works by Dorfmann and Ogden (2005); Thylander et al. (2012) and Cohen and 
deBotton (2014) among others suggested a decomposition of the stress reminiscent of the one 
obtained in Eq. (13). In some of these works it was assumed that 0” depends on F only. In 
this work, since t = T (F, E), we find that this mechanical stress component is influenced by 
the applied electric field too. We finally note that for E = 0 the well-known result of Kuhn and 
Griin (1942) is recovered. 
By substituting Eqs. (36) and (37) into Eq. (15), the polarization is 


1 oh ORT! ar& 
P=- —pdr} + . ; 39 
TdVo > ( on? i OE 0) 


where once again the summations are carried over all of the chains. The polarization Eq. (39) 
can also be viewed as the sum of two contribution. The first term is associated with the variation 
of the electrical enthalpies of the monomers as a result of the applied external electric field, and 


the second attributes to the reorientation of the chains in response to the electric field. 


4. ANALYSES OF ISOTROPIC NETWORKS WITH DIFFERENT MONOMERS 


In order to demonstrate the merit of the results developed so far, we examine a few specific 
cases, starting with the purely electrical behavior of amorphous dielectrics and tackle the electro- 


mechanical response of polymers next. To this end, we assume that the interaction between the 
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external electric field and the dielectric is solely through the dipolar molecules, and that the 
external electric field is larger than the electric fields induced on a dipole by its neighbors, thus 


deeming the internal interactions negligible. 


4.1. The local behavior of the monomer. Following Cohen and deBotton (2015), we assume 
that the dependence of the dipole vector m on the electric field can be approximated by the first 
order Taylor expansion series about E = 0. Specifically, in the monomer local coordinate system 
we assume that 


m(E) =m) + ME. (40) 


The first and second terms in Eq. (40) represent a spontaneous polarization and a linear depen- 
dence on the electric field, respectively. 

While this relation superficially appears linear, it is important to notice that it is in fact 
nonlinear. To see this, consider a single monomer, and assume a linear response m = ME, with 
M aconstant tensor. This response must be invariant under any orthogonal transformation R, i.e. 
m — RmandE — RE. This implies the trivial response y = yI. However, Eq. (40) is written in 
the monomer local coordinate system. Therefore, both mo and M are functions of é and transform 
accordingly. This enables nontrivial response and conforms to invariance under orthogonal 
transformations, at the cost of being nonlinear. This nonlinearity is conceptually analogous to 
the nonlinear strain measures that are required in elasticity for rotational invariance. It can be 
contrasted with crystalline materials in which the crystal lattice provides a “background”; there, 
one does not require invariance under any orthogonal transformation R, but only those that are 
in the symmetry point group of the crystal (e.g., James and Miiller, 1994; Marshall and Dayal, 
2014). In the crystal setting, a linear response is an adequate first-order approximation. In the 
current setting of monomers, the nonlinearity of the response is essential and is the first-order 
model. This precludes us from directly using many developed tools in statistical mechanics that 
can allow closed-form solutions for harmonic energies (e.g., Su and Purohit, 2012). 

The corresponding electrical enthalpy of a dipole oriented along é is (Blythe and Bloor, 2008; 
Cohen and deBotton, 2014) 

h=-m-E. (41) 


For later reference we note that h # h(F). Next, we list the three specific models for the local 


relation (40) that will be utilized in the forthcoming analyses. 
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The first corresponds to a spontaneous dipole, or a rigid dipole with a constant magnitude & 


(Davidson 1962; Blythe and Bloor 2008), 


ms =Sé. (42) 


Following a concept proposed by Stockmayer (1967), Cohen and deBotton (2015) introduced the 
second model of a uniaxial dipole whose magnitude depends on the electric field. Specifically, 
2 

my = U8 OFF, (43) 
where we formally normalize by k T so that Ry has the units of a dipole. 
Another type of dipole discussed by Stockmayer (1967) was represented by Cohen and 


deBotton (2015) as the transversely isotropic (TI) model. In this case the dipole 


mr; = 1 871 (1-202), (44) 
2k 
is perpendicular to € and 87 has the units of a dipole. 

For convenience, throughout this work we set Ry = Kr; = K. This will ensure that three 
dielectrics composed of a random and uniform distribution of spontaneous, uniaxial and TI 
dipoles admit the same behavior in the limit of infinitesimal deformations and small electric 


fields. 


4.2. The purely electrical response of amorphous dielectrics. Consider an amorphous di- 
electric with no = am dipoles per unit referential volume subjected to an electric field E. We are 
concerned with the purely electrical response, and to eliminate the mechanical constraint (20) 


we set tT = 0 in Eq. (39). The resulting polarization of the dielectric is 
P=n0 [ mpar. (45) 


where Eqs. (30) and (31) result in the Boltzmann distribution p = 5 exp (22) with Z = 
J exp (2F) dl and J = 1. Next, we consider the three local behaviors of the dipoles. 
With Eq. (45), the polarization of the dielectric is 


Ps =m KL (wv) E, (46) 
where w = ae and 2 (w) = cothw — + is the Langevin function. Note that as E — oo (or 


w — oo) all dipoles tend to align with the electric field and therefore Ps — no& E. At small 
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electric fields the Taylor expansion series up to first order in w is 


8 
Ps~ rw (47) 


where at E = 0 (or w = O) we find that Ps = 0 due to the initial random distribution of the 


dipoles. The corresponding expression for the initial susceptibility 


. noK2 
~ 3kT e0 


X0 (48) 


is in agreement with the results of Debye (1929), Frohlich (1949), Davidson (1962) and Blythe 
and Bloor (2008). 

We consider next a material with a random distribution of uniaxial dipoles. From Eq. (30), 
the partition function is 


(49) 


where D (w) = exp (-w?) i exp (:?) dt is the Dawson function or the Dawson integral. There- 
fore, 
exp (—w” sin? (6;)) . (50) 


= W 
PU An Do) 


Integration according to Eq. (45) yields 


K 
Py = = nu ob, (61) 
where ny = aoe . In the limit of small electric fields ny = = and P = mS WE. At the other 


limit, as E — oo (or w > ©), ny — 2 and all the dipoles tend to align with the electric field. 
Consequently, Py — no x E. We note that 7y (w) is a smoothly varying bounded function 
between 5 at w =O and2 asw > o~. 
Next we consider the case of a TI dipole Eq. (44). From Eq. (30), the partition function is 
(2 n)>/? w W 
ZrI = —| Erf |—], 52 
TI pms ae ee Va (52) 


where Erf (x) is the well-known error function, and the corresponding PDF is 


= = exp ( 
3/2 wo 
(2.2)3/2 Erf ( +) 


By means of Eq. (45) the resulting polarization is 


2, 
Pri _ cos? co] (53) 


no ® 
Pr) = > nr wk, (54) 
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Ficure 1. The normalized polarization versus w. The continuous, dotted, dashed 
and dot-dashed curves correspond to the polarization according to the linear, the 
spontaneous, the uniaxial and the TI dipole models, respectively. 


exp|-%- 
where nr; = 1 - 4+ + a2 cal) Again we note that at small electric fields n77 = 4 and 


wErt(4) 


Pr; = mS wot, and as E > o (or w > ©9), nr; — 1 and Pr; > + no SE. The function 
nT1 (w) is a smoothly varying bounded function between w = 0 and w > ov. 

The variations of the polarization with the electric field according to the three models together 
with a linear model with a constant susceptibility are compared in Fig. (1), where the polarization 
in the direction of the electric field normalized by no is plotted versus the dimensionless 


quantity w. The continuous curve corresponds to the commonly assumed linear model 


Pe= me i (55) 


This model does not take into account the evolution of the dipoles configuration in response to 
the external field. At small fields all three models agree with the linear one. The dotted curve 
represents the predictions according to the spontaneous model and at large electric fields, as the 
dipoles orient along the field, the overall polarization saturates. The dotted and the dot-dashed 
curves depict the predicted polarizations according to the uniaxial and the TI dipole models, 
respectively. We note that these curves admit a form of bilinear curves with distinct initial and 
final slopes. Due to the choice of the model constants, as w — oo the ratio between the two 
predictions approaches 2 (in accordance with the limiting values of ny and 777). 

We conclude this section noting that by taking a linear combination of the three models 
described in section (4.1), dipoles of various orientations and intensities with respect to the 


electric field and the direction of the molecule é can be considered. 
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4.3. The electro-mechanical response of polymer networks. Consider a polymer with No 
chains per unit referential volume where each chain is composed of n dipolar monomers of 
length /. At the reference configuration the vector that connects the two ends of the k-th 
chain is aligned with the direction t” and has a length Yn/ (Treloar, 1975). At the reference 


configuration the chains are randomly oriented and uniformly distributed such that 


(fo) = 0, (56) 
and 
1 
(Fo @ Fo) = 3 I, (57) 
where 
1 yee) 
= —__ dl = : 58 
(*) No dV [oe No dVo oy) 


denotes the average of a chain related quantity e. Throughout this work, the averages are 
determined by application of the numerical micro-sphere technique that was employed, for 
example, by Miehe et al. (2004), Thylander et al. (2012), and Cohen et al. (2015). Relying on 
the experience gained in those works, we consider 42 integration points. 

In accord with the works of Treloar (1975) and Flory (1953), we assume that when subjected 
to an electric field E and a deformation gradient F the chains undergo affine deformations, i.e. 
for the k-th chain 

r) = val FR”, (59) 


We note that, while from this relation it follows that the electric field does not affect directly the 
distribution of the chains (i.e., r“ 4 r“ (E)), at the monomer level the electric field can lead 
to rearrangement and redirection of the dipolar monomers. 


Next, we define the dipole of the k-th chain 
m\*) = yim” =n [mp dI, (60) 
i=l 


where the summation is carried over all the dipoles composing the chain and, assuming a large 
n, we execute the summation in terms of an appropriate integral. According to Eq. (39) the 


polarization is 
N 
P = — (mc). (61) 
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Since h # h(F), the mechanical stress Eq. (38) reduces to 


kT 
om — No vn 


5 (t @fo) F’. (62) 


Utilizing Eq. (61), we first determine the polarization of a polymer whose chains are made of 


spontaneous dipoles Eq. (42), and find that 
Non a N K 
p=—— ( | par) =" ray =0, (63) 


where Eqs. (60), (32), (59) and (56) are used. This result can be explained as follows: since 


the magnitude of the spontaneous dipoles is fixed, the electric field can only influence their 
local orientations. Due to the mechanical constraint Eq. (20), the directions of the end-to-end 
vector and the dipole of any polymer chain coincide. Upon integration from the chain to the 
macroscopic level, we follow the affine deformations assumption and find that the sum of the 
dipoles of any pair of chains lying along the same vector but pointing in opposite directions 
vanish. 

The mechanical component of the stress can be derived from Eq. (13). We note that while the 
individual dipolar monomers would energetically favor orienting themselves along the electric 
field, they are constrained from doing so by the fact that the end-to-end vector is constrained 
by F. While this constraint does not allow a net dipole moment to develop, it also leads to a 
mechanical stress that physically reflects the inability of the monomers to be oriented favorably. 
The emergence of mechanical stress in response to the affine deformation constraint will be 
discussed further in the following subsection. 

We are interested in the electro-strictive properties of polymer networks, and therefore in the 
following we consider chains made out of dipolar monomers whose magnitudes do depend on 
the electric field. Particularly, we examine networks with uniaxial Eq. (43) and TI Eq. (44) 
dipoles. Accordingly, consider a polymer composed of long chains of uniaxial dipoles. The 


PDF Eq. (31) of the distribution of the dipolar monomers in a chain is 


p (8h) = 5 exp(t-E +0 (2-£)’). (64) 
Due to the assumed large number of dipolar monomers in a chain there is a finite range of 


deformations for which 


r or 1 
gc a ee a Sa 1 
aoa we fo - Cfo « 1, (65) 
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where Eq. (59) is used and C = F’ Fis the right Cauchy-Green deformation tensor. Consequently, 
by approximating the PDF in Eq. (64) with its second order Taylor expansion and employing 
Eq. (32) we derive the approximation 


4 
2—Nu 


1 
Ty = (I ~ayke E) e Fio, (66) 


where ay = Si A closed-form estimate for the chain’s dipole can be obtained from Eq. (60) 


nKw 1 1 2 
m\*? = 5} ug (1+ 75 (5 (Cau - bu) tu TU + (3by —S5ay) (tu -E) )1+ au ty @ tv) f. 
(67) 
where by = w (1 - Ty), To calculate the polarization we make use of Eqs. (66), (56), (57) and 


(59) and arrive at the expressions 


ae ae a 
(ry ty) = (—_ (tr (b) + (az — 2a) E- bE), (68) 
" 1/ 4 \/ 5 
((ru-8) ie 3h (—_] (1 - ay) E - bE, (69) 
and 
1( 4 \ 
(ty @ Tu) = (—_ (b-ay (E®bE+bE@E)+a,£-bEE@E), — (70) 
3n 2-nu 


where b = FF’. Subsequently, Eq. (61) becomes 


2 
py = eRe rs (—_] (5 ((au = bu) 1) + C5 B- bE) T+ ay b- C2] f, 


i 2 w*n \2—nu 
(71) 
where 
Ci? = (ay - bu) « (ay - 2av) + Bbu - Sav) (1- au)’, (72) 
and 
C\? = aj, (R@ bE +bE @£) +a}, (E- bE) Lek. (73) 


We point out that a Taylor series expansion around w — 0 for the two functions reveals that 
c. ~wt+O (w°) and CO ~wi+O (w°), whereas the coefficients of tr(b) and b are 
proportional to ~ w* + O (w*). Therefore, for w < 0.3 the polarization may be approximated 
via 


7 No &n 


a (I~ 5 crt 3b)) of. (74) 
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The zero order term in Eq. (71) is isotropic and the associated contribution to the polarization 
field is aligned with the external electric field. This term corresponds to the initial macroscopic 
susceptibility of the polymer and is identical to the result obtained in Eq. (48) for amorphous 
dielectrics with no = Non. The first order term, which is proportional to ‘, stems from the finite 
length of the chain and has an influence on the polarization as the polymer deforms. We compare 
this result to the one obtained in the work of Cohen and deBotton (2015), where a deterministic 
analysis in which the entropy was not accounted for was introduced. While the zero order term 
from both analyses is identical, the first order term in the earlier work is different from the one 


obtained herein and proportional to This discrepancy may be explained as follows: even 


though the analysis in Cohen and es (2015) is based on valid chain configurations, these 
are not the most probable ones. We argue that the current analysis is more physically sound since 
the response of polymers is entropy-driven and the present work considers the configurations 
that maximize the entropy. 


The stress is computed via Eq. (13), where the mechanical stress Eq. (62) is calculated with 


Eq. (66), 


No kT 4 
Ou = 


a7 5 (I-wEe8))b+E@Py +0", (75) 


where Eqs. (57) and (66) are employed. By setting E — 0 (or w — 0), the above expression 


reduces to the standard neo-Hookean stress, where 
= NokT, (76) 


is the shear modulus (Treloar, 1975). 
By following similar steps for a polymer with chains that are composed of TI dipoles Eq. (44), 
we approximate tT for long chains as 


2 2 - 3nr1 
= —(1+——“"_£ a k| —F 77 
— | "2 @n-D . Fc we eg 


The expression for the chain dipole in this case is similar to the one in Eq. (67), where ay and by 
are replaced with az; = -_ —3 and br; = w* (nz; — 1), respectively. The resulting expression 


for the polarization is 


& 1 a ago 
pry = not (14 (—] (5 (er = orn tr(b) + Cj B- BE) I+ arb - CF] E. 


J 2 
(78) 
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Here the expressions for Co and C ce are given in Eqs. (72) and (73), respectively, and where 
here too ay and by are replaced with ary; and br,, respectively. Following arguments similar to 
the ones led to Eq. (74), we argue that for w < 0.3 the polarization may be approximated via the 


Taylor expansion series for small w 


No &n 1 
Pry © I tr (b) I- E. 79 
TI ( + T9, (er ®) 3)) (79) 
The total stress is 
NokT [{ 2 
OTI= — (I-ar;E@k) b+E@P7,+0", (80) 
3J)— \nri 


where Pv; is given in Eq. (78). 


5. APPLICATIONS TO POLYMER FILMS 


The polarization and the stress in Eqs. (61) and (13) depend on three key microscopic factors: 
the type of the dipolar monomers composing the chains (i.e., uniaxial or TI dipoles), the 
local dipole constant (8, or X77), and n the number of dipolar monomers in a chain. These 
microscopic factors will dictate the macroscopic response and are therefore prominent for the 
electromechanical coupling at the macroscopic level. We investigate the impact of these three 
microscopic factors on the macroscopic behavior. To this end, we examine thin layers of 
polymers with different dipole types, dipole constants and chain lengths that are subjected to the 


equi-biaxial isochoric deformation 
1 
F=1(I-Rek)+ Fok. (81) 


This type of loading is common in experimental settings where the faces of the film are covered 
with flexible electrodes of negligible stiffness that are charged with opposite charges such that 
the electric potential difference (AV) induces an electric field E = EF E across the film (see 
Fig. (2)). Throughout, we assume that the shear modulus of the polymer in its initial unloaded 


configuration is = 10° Pa, and deduce No, the density of the chains, from Eq. (76). 


5.1. The role of the local dipole constant. In order to examine the effect of the constant 8 we 
compare the behaviors of two polymers with initial susceptibilities Ne = 3.7 and 5 =o, 
The chains composing these polymers contain n = 100 uniaxial monomers. We choose the value 
of ye to be equal to the electric susceptibility of the commercially available VHB 4910, which 


is commonly used in electro-mechanical actuation experiments (Kofod et al. 2003; McKay et 
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Undeformed configuration Deformed configuration 


* flexible electrodes 
BE 


Ficure 2. A sketch of a polymer film subjected to an electric field and a me- 
chanical deformation. 


al. 2009; Di a et al. 2012; Qiang et al. 2012). As before, the dipole constants are 8\ = al 
and he = a a where §; and > are determined from Eq. (48) with no = Non. 

Fig. (3) shows the distributions of the dipolar monomers belonging to a single chain whose 
end-to-end vector lies on the plane spanned by FE and Y. In these figures the length of the radius 
vector to a point on the curve is proportional to the number of dipolar monomers aligned with this 
vector. We note that an isotropic distribution of monomers in a chain whose end-to-end vector 
vanishes is represented by a circle with a radius 7-. For convenience, in the polar distribution 
plots we normalize n by this isotropic distribution such that if the directional distribution of 
the monomers in a chain is isotropic, the radius of the corresponding circle will be unity. To 
gain some quantitative information from these plots we also draw light circles with normalized 
radii 1/2, 1, 3/2 and so forth. Specifically, the number of monomers aligned with a ray from 
the origin to an intersection of these circles with the curves for the directional distribution is 
equal the normalized radius of the circle multiplied by the factor 7-. Moreover, the amount by 
which the directional distribution is deviating from the unit circle provides an idea regarding its 
anisotropy. 

Figures (3a) and (3b) correspond, respectively, to the two polymers with low and high dipole 
constants, where both polymers are subjected to a transverse stretch ratio 2 = 5. In both 
figures the chain’s end-to-end vector is aligned with the direction of the electric field, and the 
continuous, dashed and dot-dashed curves correspond to different intensities of the applied 
field. In the absence of an electric field and under the imposed stretch, the affine deformation 
assumption results in a significant shortening of the end-to-end vector of these chains and thus 
the distribution of the monomers is nearly uniform. When an electric field is imposed, since the 
response of a uniaxial dipole is quadratic in é, it causes the monomers to rotate toward either 
E or -E. Indeed, in Fig. (3) we note the tendency of the monomers to align with the electric 


field, and the higher the electric field is the stronger this tendency becomes. Moreover, we 
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(a) ~ (b) 


Ficure 3. The distribution of the uniaxial dipolar monomers in chains aligned 
with the electric field belonging to polymers with (a) yo = 3.7 and (b) yo = 37 
subjected to a biaxial stretch with A = 5. 


EB E = 100MN 


Ficure 4. The distribution of the uniaxial dipolar monomers in chains whose 
end-to-end vectors are perpendicular to the electric field belonging to polymers 
with yo = 37 subjected to a biaxial stretch with A = 5. 


also note that the rate at which the monomers align with the field strongly depends on the local 
dipole constant, implying that at the chain level as well as at the macroscopic level the overall 
polarization is not linear with &. 

Fig. (4) is similar to Fig. (3) but depicts the polar distribution of the dipolar monomers 
belonging to a chain whose end-to-end vector lies along the Y-direction. Since the variations in 
the distributions of the monomers in the polymer with yo = 3.7 are small, here and subsequently 
we show the polar distributions of the monomers polymers with yo = 37 only. Due to the 
imposed deformation, this chain is stretched and the monomers are inclined with the principle 
stretch direction. The electric field tends to reorient the monomers perpendicularly to increase 
their average projection in the field direction. Once again, the larger the field is, the stronger the 


monomers tendency to reorient is. 
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Regardless of the chain’s end-to-end direction, we note in Figs. (3) and (4) the tendency of the 
monomers to align with the electric field. This, as demonstrated in Fig. (Sa), results in an increase 
of the macroscopic polarization with the electric field. In Fig. (5a) the normalized component 
of the polarization vector P = Fe, with Po = yoeoEH, is plotted against E for polymers with 
the two susceptibilities. Accordingly, at E = 0 the limit P = limgo re is taken. Note 
that due to the dependence of the coupled properties on the deformation, even in the absence 
of an electric field P # Po. Fig. (5a) illustrates that the macroscopic susceptibility depends 
on the electric field and the intensity of this dependence is amplified in polymers with higher 
microscopic dipole constants. This can be explained with the aid of the PDF for the distribution 
of the monomers as follows: the dipole constant has a strong influence on the rearrangement 
of the dipolar monomers in response to the electric field. Thus, on top of the increase in the 
overall polarization due to the local behavior of the monomers, the freedom the monomers have 
to rotate further enhances this tendency. In accordance with this observation we note that even 
though the ratio between the initial susceptibilities of the two polymers shown in Fig. (5a) is 
10, at high electric fields the ratio between the polarizations of the two polymers is much larger. 
For example, at E = 100 MV we find that ee = 1.2. The foregoing analysis further 
contradicts the common assumption that the susceptibility of a polymer is constant (Dorfmann 
and Ogden, 2005; Shmuel and deBotton, 2013). In passing, we note that the second term in 
Eq. (75), namely the polarization stress E @ P, exhibits the same dependence on the electric field 
as the one discussed in connection with Fig. (5a). 

We recall that the mechanical stress 0” Eq. (62) is one of the three components of the 
stress in Eq. (13). Fig. (5b) depicts the normalized stress a, = a aR 


where N is a unit vector and a) =a" (A, E = 0) is the purely mechanical stress (Kuhn and 


as a function of E, 


Griin, 1942; Treloar, 1975). Following the curves for this stress component in Fig. (5b) we 
observe that in polymers composed of uniaxial dipoles of decreases with the electric field 
while oy increases. To explain the observed variations of this stress component we recall 
that, as discussed in connection with Figs. (3) and (4), the dipolar monomers tend to reorient 
themselves along the electric field. From a mechanical-geometrical viewpoint this rotation will 
result in lengthening the projections of the monomers along the electric field and shortening 
their projections in the transverse direction. This local response prompts similar tendencies of 
the corresponding projections of the chains end-to-end vectors. However, since according to 


the affine deformation assumption the deformations of the end-to-end vectors are constraint, 
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Ficure 5. (a) P and (b) a7, and oy versus E for polymers with uniaxial dipolar 


monomers and two dipole constants subjected to a stretch A = 5. 


the inability of the monomers to freely reorient manifest in terms of a compressive stress 
variation in the direction of the electric field and a tensile stress variation in the transverse 
direction. We note that while the applied electric field triggers these variations in the stress, this 
stress component stems from a mechanical constraint and not directly from the electric forces 
applied on the dipolar monomers. We also note that for the uniaxial dipolar monomers the 
variation of this stress component with the electric field opposes the corresponding variation of 
the polarization stress. A comparison between the two pairs of curves for the polymers with 
different dielectric constants demonstrates that the dependence of these stress variations on the 
local dipole constant is non-linear, and it is more pronounced in the polymer with the larger 
microscopic dipole constant. 

We examine next the effect of the mechanical loading on the evolution of the micro-structure 
in chains with 100 uniaxial dipolar monomers. We consider a polymer with an electric sus- 
ceptibility yo = 37 subjected to an equi-biaxial deformation. Figs. (6a) and (6b) illustrate 
the distribution of the dipolar monomers in a chain whose end-to-end vector lies along the 
E-direction for A = 1 and A = 5, respectively. The dotted and the continuous curves represent 
the polar distributions of the monomers in the absence of an electric field and when E = 50 MV 
respectively. Since the initial length of the end-to-end vector is fairly large in comparison with 
the chain contour length r = 0.1 n/, the referential distribution of the monomers is asymmetric 
with respect to the plane transverse to the end-to-end vector. Application of electric field causes 
rotation of the monomers toward E (or —E) but, since the length of the end-to-end vector is fixed, 
the distribution is still uneven with respect to the transverse plane. As the polymer deforms, 


the affine deformation assumption imposes a shortening of the end-to-end vector and hence a 
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Ficure 6. The distribution of the uniaxial dipolar monomers for chains aligned 
with the electric field belonging to polymers with yo = 37 subjected to a biaxial 
stretch (a) A = 1 and (b) A = 5. The dotted and the continuous curves correspond 
toE=0 MV and E = 50 wv respectively. 


reduction in the ratio >. As a result a more symmetric distribution of the dipolar monomers is 
exhibited in Fig. (6b). We point out that the decrease in the ratio + increases the freedom of the 
monomers to reorient, and thus enhances the influence of the electric field on the distribution. 

Fig. (7) depicts the variations in the distribution of monomers of a chain whose end-to- 
end vector lies along the Y-direction. In this case, at the referential state the distribution is 
symmetric about the Y-axis. Due to the quadratic dependence of the uniaxial dipoles on é and 
the symmetric distribution with respect to the Y-axis, the rotation of the monomers toward E 
and —E is symmetric and hence the distribution remains symmetric when the chain is subjected 
to electric field. Figs. (7a) and (7b) show the distributions of the uniaxial dipolar monomers 
belonging to a polymer subjected to E = 50 MV and biaxial stretches with A = 1 and A = 5, 
respectively. We inspect the ratio between the maximum radii with (continuous curve) and 
without (dotted curve) an applied electric field, and note that this ratio decreases as the polymer 
is stretched. We therefore conclude that the influence of the electric field on the rearrangement 
of the dipolar monomers in a chain diminishes with the stretch. We stress that in general, under a 
volume preserving deformation, there are chains whose monomers gain degrees of freedom and 
become more sensitive to the application of E and chains in which the mechanical constraints 
become stricter, limiting the rotation of the monomers, and thus diminishing the effect of the 
electric field. 

Next, we examine the effects of the deformation on the macroscopic response of polymers 
with different electric moduli. Fig. (8a) depicts the dependence of P on the planar stretch A? of 


two polymers with different susceptibilities. As the polymer deforms, the mechanical constraints 
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Ficure 7. The distribution of the uniaxial dipolar monomers in chains whose 
end-to-end vectors are perpendicular to the electric field belonging to polymers 
with yo = 37 subjected to a biaxial stretch (a) A = 1 and (b) A = 5. The dotted 
and the continuous curves correspond to E = 0 and E = 50 MV respectively. 


on chains with end-to-end vectors on the Y — Z plane lead to a rotation of the uniaxial dipolar 
monomers away from FE and, as a result, the overall polarization decreases. We note that the 
curves describing the predicted polarization of the two polymers are roughly parallel. This is 
anticipated since the reorientation is due to a mechanical deformation at a fixed electric field 
and hence the influence of the susceptibility is fairly small. 

Fig. (8b) shows the variation of the mechanically-dominated normalized stress components 
&" along the E and Y directions versus A* for the same two polymers. Before we proceed, we 
stress that o”” strongly depends on the deformation and, as expected, admits the lock-up effect 
due to the limiting chain extensibility of the polymers. However, since in this work we explore 
the coupled behavior of the polymer, we examine the difference in the dependence of the stress 
on the deformation with and without the application of an electric field. For this reason we show 
the variations of the normalized quantity &-””. Note that for the two polymers the two components 
of the normalized stress @” are roughly constant. This implies that under a prescribed electric 
field, the trend of the mechanically-dominated stress is almost identical to the one obtained in 
the absence of the electric field. Nonetheless, due to the effect of the electric field at A? = 1, 
oO” #1. 

Next, we examine two polymers with transversely isotropic monomers. Their initial suscep- 
tibilities are 3.7 and 37 as before, and the chains are made up of 100 dipolar monomers. To 
highlight the response of these monomers to the electric field we examine first two specific 
chains of the polymer with yo = 37. Ina manner similar to Figs. (3) and (4), Figs. (9a) and (9b) 


present polar plots of the distributions on the plane spanned by E and Y of the dipolar monomers 
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Ficure 8. (a) P and (b) oO, and oy versus the planar stretch A* for chains 
composed of uniaxial dipolar monomers with two dipole constants subjected to 
an electric field E = 50 MY 
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Ficure 9. The distribution of the TI dipolar monomers in chains whose end-to- 
end vectors are along (a) the E and (b) the Y directions belonging to a polymer 
with yo = 37 under a stretch with A = 5. 


in chains whose end-to-end vectors are along the E and the Y directions, respectively. These 
plots show that the electric field causes the monomers in both chains to rotate away from the 
direction of the electric field. This is a consequence of the local behavior of the TI dipole. 
Thus, since the dipoles are perpendicular to é, the tendency of the electric field to reorient them 
in order to increase the projection of their electric dipoles in the direction of E results in the 
observed behavior. 

Fig. (10a) depicts the normalized polarization P as a function of E. The legend of this figure 
is identical to the one used in Fig. (5a). The polarization increases non-linearly with the electric 
field. In contrast with polymers with uniaxial dipoles, the pre-stretch results in redistribution 
of the monomers such that their dipoles rotate toward the electric field, resulting in an initial 


enhancement of the polarization. Fig. (10b) shows the mechanical stress oy, as a function of 
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Ficure 10. (a) P and (b) ao, and ay versus EF for polymers with TI dipolar 


monomers and two dipole constants subjected to a stretch A = 5. 


E, where we heed that the trends of the variations in 0 and oy’ are opposite to the ones 
revealed for polymers with uniaxial monomers. Thus, due to the shortening of the monomer 
projection in the E-direction, the variation in the mechanical stress along the direction of the 
electric field becomes more tensile, while in the perpendicular direction the stress becomes 
more compressive. We note that the dependence of the mechanical stress variations on E in 
polymers with TI monomers result in enhancement of the polarization stress in the direction of 
the electric field. This implies that polymers with TI dipoles may be preferable in applications 
for EAP actuators. 

Figs. (1 1a) and (1 1b) depict the polar distributions of the 100 TI dipolar monomers composing 
a chain whose end-to-end vector is along the E-direction on the plane spanned by E and Y for 
various stretch values. Analogous plots for a chain whose end-to-end vector is along the Y- 
direction is shown in Figs. (12a) and (12b). When excited by an electric field (continuous curves) 
the TI dipolar monomers of both chains rotate away from the direction of the electric field in 
order to orient the dipoles along its direction. At high stretch ratios a chain with an end-to-end 
vector in the direction of the electric field shortens, resulting in a more symmetric distribution 
of the monomers orientations. This also provides the monomers more freedom to reoriet with 
the electric field. At the other hand, the end-to-end vectors that are perpendicular to the electric 
field elongate at A = 5 and hence the monomers ability to reorient with the field diminishes. 
This can be seen from the proximity of the continuous and the dashed curves in Fig. (12b). 

In Figs. (13a) and (13b) P and the components of 0” are plotted as functions of A? for polymers 
composed of 100 TI dipolar monomers with two different susceptibilities, respectively. As 


discussed previously with relation to Figs. (12a) and (12b), as the planar deformation increases, 
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Ficure 11. The distribution of the TI dipolar monomers in chains aligned with 
the electric field belonging to polymers with yo = 37 subjected to a biaxial pre- 
stretch (a) 2 = 1 and (b) A = 5. The dotted and the continuous curves correspond 
toE=0 MV and E = 50 wv respectively. 


(a) (b) 


Ficure 12. The distribution of the TI dipolar monomers in chains whose end- 
to-end vector is perpendicular to the electric field belonging to polymers with 
Xo = 37 subjected to a biaxial pre-stretch (a) A = | and (b) A = 5. The dotted 
and the continuous curves correspond to E = 0 MV and E = 50 MV respectively. 


the freedom of the monomers belonging to chains whose end-to-end vectors lie on the plane 
transverse to the electric field to rotate becomes limited. However, thanks to the local behavior 
of the TI dipoles, the planar deformation causes these dipoles to align with the electric field 
and, consequently, to increase the polarization. As expected, this effect is more pronounced in 
polymers with higher susceptibilities. 

The trend of the mechanically-dominated stress components o7; and oy; shown in Fig. (13b) is 
similar to the one predicted for polymers composed of uniaxial dipoles (Fig. (8b)). However, as 
mentioned in the context of Fig. (10b), the electric field increases the mechanically-dominated 
stress component in the direction of the electric field and decreases the component in the 


transverse direction. 
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Ficure 13. (a) P and (b) oO, and ay’ versus A? for polymers with TI dipolar 
monomers and two dipole constants subjected to an electric field E = 50 MV 


5.2. The role of the length of the chain. The influence of the number of dipolar monomers 
per chain on the overall response is examined next. We consider four polymers with chains 
composed of 25, 100, 500 and 1000 uniaxial monomers and an initial susceptibility yo = 3.7. 
Figs. (14a) and (14b) depict P as a function of the electric field in polymers subjected to biaxial 
pre-stretch with A = 3 and A = 5, respectively. The dashed and continuous curves correspond 
to the long-chain approximation (71) and the associated full model predictions, respectively. 
The polarization of polymers with n = 500 and n = 1000 is almost constant with the electric 
field under the imposed stretches. For these polymers the long-chains approximation provides 
excellent predictions. For polymers with chains composed of 100 monomers at A = 3 this 
approximation still holds. If the polymer is further stretched up to 2 = 5 the long-chain 
approximation slightly overestimates the polarization. 

In general, Figs. (14a) and (14b) demonstrate that the polarization increases with the electric 
field. This effect is more pronounced in short-chain polymers as evident from the results for 
n = 25. This short-chain polymer locks-up at 2 = 5, and hence up to E = 100 has a very low 
polarization (P in the range between 0.4 and 0.44), which is not shown in Fig (14b). 

Another important observation that can be drawn from these plots is that the polarization of 
polymers with uniaxial dipoles increases with the length of the chain. In order to explain this 
observation, consider two chains with 25 and 1000 monomers whose end-to-end vector lies on 
the plane perpendicular to E. In the limit A — 5 the short chain is fully stretched, and its 
uniaxial dipoles are oriented normal to the electric field. Consequently the dipole of this chain 
is practically zero. Conversely, at this stretch ratio the monomers of the long chain still possess 


freedom to reorient and some have projections along the electric field, rendering a non-zero 
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Ficure 14. P versus E according to the proposed model for polymers with 
Xo = 3.7 and chains containing different numbers of uniaxial dipoles subjected 
to biaxial stretch loading with (a) A = 3 and (b) 2 =5. 


chain dipole. These two extreme examples demonstrate the higher sensitivity of polymers with 
shorter chains to the applied stretch. 

It was shown in connection with Fig. (14) that the macroscopic behavior of polymers strongly 
depends on 1 when the deformation is held fixed and the electric field is varied. Next, we 
carry the opposite analysis and fix the electric field while varying the deformation. Figs. (15a) 
and (15b) depict P as a function of the planar stretch with E = 1 MV and E = 50 MV for four 
different values of n, respectively. The polymers are made out of uniaxial dipolar monomers and 
their initial susceptibility is yo = 3.7. We once again employ the long-chains approximation 
Eq. (74) for the cases n = 1000 and n = 500. Due to the normalization of P - E, the curves 
obtained for E = 1 MV and E = 50 MV for the polymers with the long-chain approximation are 
practically identical. The long-chains approximation holds nicely even for polymers with 100 
monomer chains. Furthermore, it is again demonstrated that the variation of the polarization is 
proportional to the chain length. 

Figs. (16a) and (16b) show P as a function of E for four polymers with chains of different 
lengths composed of TI dipolar monomers subjected to biaxial stretch with 2 = 3 and A = 5, 
respectively. Here, the stretch results in aligning of the dipoles of the chains whose end-to-end 
vectors are perpendicular to E and consequently two main effects are observed. The first is an 
increase in the polarization with the stretch and the second is that short-chain polymers exhibit 
higher growth rate of the polarization with the stretch. 

Lastly, we examine the dependence of the responses of polymers with different chain-lengths 


composed of TI dipolar monomers to the applied deformation. The predicted normalized 
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Ficure 15. P versus A? for polymers with yo = 3.7 and chains composed of 
uniaxial dipoles with different chain lengths subjected to (a) E = 1 MV and (b) 
E=50™. 
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Ficure 16. P versus E according to the proposed model for polymers with 
Xo = 3.7 and chains containing different numbers of TI dipoles subjected to 
biaxial stretch loading with (a) A = 3 and (b) A =5. 


polarization P as a function of the planar stretch for applied electric fields E = 1 MV and 
E = 50 MV are shown in Figs. (17a) and (17b), respectively. The predicted polarization 
increases with the stretch, as explained in connection with Fig. (13a). In this case too the 


long-chains approximation holds for chains with 100 TI dipolar monomers. 


6. CONCLUSIONS 


This work introduces an entropy-based multiscale analysis of the electromechanical coupling 
in dielectric elastomers. Initially, the variational principles for the electrical enthalpy-density 
and the entropy-density are written in accordance with the first and the second laws of thermo- 
dynamics. Next, the conformations that a polymer chain can occupy are considered. Following 


common practice, we assume that the configuration of the chain that maximizes the entropy is 
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Ficure 17. P versus A? for polymers with yo = 3.7 and chains that are composed 
of TI dipoles with different chain lengths subjected to (a) E = 1 MV and (b) 
E=50™. 


the one corresponding to the most probable distribution of the dipolar monomers. Next, the 
general expressions for the polarization and the stress that develop in the dielectric are derived. 

To determined the macroscopic coupled response of specific polymer networks three models 
for the local behaviors at the monomer level are considered: a spontaneous dipole (Blythe 
and Bloor, 2008), a uniaxial dipole and a TI dipole (Stockmayer, 1967; Cohen and deBotton, 
2015). We first consider the purely electrical response of amorphous dielectrics. The predicted 
polarization of the dielectric in response to an electric field according to the three models is 
determined and compared with a macroscopic approach that assumes a constant permittivity. 
We find that none of these models agrees with the assumed linear behavior, and while the 
spontaneous dipole model leads to saturation of the polarization, the other two models admit a 
bi-linear response at the macroscopic level. 

Next, a long-chain approximation for the polarization and the stress in polymers with chains 
composed of a large number of dipolar monomers is determined. We reveal that the macroscopic 
coupled response of the polymers crucially depends on three key parameters — the type of 
the dipolar monomers composing the chain, the magnitude of the microscopic dipoles and 
the number of dipolar monomers in a chain. Accordingly, a thorough investigation of the 
influence of these three parameters is carried out by tracking the variations of the directional 
distributions of the monomers along the chains. In passing, by comparing with the full model, 
we also demonstrate that the closed form long-chains approximation is valid for a wide range of 
polymeric chains and stretch ratios as large as 3. A particularly interesting result concerns the 


dependence of the mechanical stress on the electric field. We find that this dependence, which is 
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commonly neglected in phenomenological models, can have a large impact on the macroscopic 
response. Moreover, these electric field induced variations in the mechanical stress may enhance 
or oppose the polarization stress stemming from the interaction of the dipolar monomers with 
the field. 

The findings of this work shed light on the relations between the molecular structure of 
polymers and the origin of their coupled electromechanical behavior. In turn, the analysis 
presented herein can pave the path to the enhancement of the poor electromechanical coupling 
in dielectric elastomers that hinders their widespread applications. Finally, while beyond the 
scope of this paper, one can also derive an electroelastic energy density starting either from 
the partition function or through transformation of the macroscopic thermodynamic quantities 
H and S. Such an energy density would provide microscopically-motivated constitutive input 
to calculations such as ones carried out by Xie et al. (2008); Rudykh et al. (2013); Tian et al. 
(2012); Galipeau and Ponte Castaneda (2012); Gei et al. (2013); Spinelli et al. (2015), among 


others. 
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Abstract 


Nematic liquid crystals composed of rod-like molecules have an orientational elasticity that ac- 
counts for the energetics of the molecular orientation. This elasticity can be described by a unit vector 
field; the unit vector constraint interacts with even fairly simple boundary conditions to cause discli- 
nation defects. Disclinations are entirely a topological consequence of the kinematic constraint, and 
occur irrespective of the particular energetic model. Because disclinations are topological defects, 
they cannot be regularized by adding higher gradients, as in phase-field models of interface defects. 
On the contrary, the higher gradient terms would cause even greater singularities in the energy. In this 
paper, we formulate an integral-based nonlocal regularized energy for nematic liquid crystals. Our 
model penalizes disclination cores and thereby enforces a finite width, while the integral regulariza- 
tion ensures that the defect core energy is bounded and finite. The regularization at the same time 
tends to the standard gradient-based energies away from the disclination, as well as building in the 
head-tail symmetry. We characterize the formulation in its ability to describe disclinations of various 
strengths, and then apply it to examine: (1) the stability and decomposition of various disclinations, 
and the competition between bend and splay energies in determining the relative stability of inte- 
ger and half-integer disclinations (2) the coalescence of a +5 and —4 disclination pair; we find the 
disclinations do not move at the same velocities towards each other, suggesting that the asymmetry 
of the director field plays a dominant role despite the equal-and-opposite topological strengths of the 
disclinations. 


1 Introduction 


Nematic liquid crystals are composed of rod-like molecules that have no positional ordering but tend 
to align along the same direction. There is therefore an energetic penalty — an orientational elasticity 
— for spatial variations in the direction. In a continuum mesoscale description, the kinematics of the 
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orientational ordering can be described by a unit vector field n— the director field — that tracks the average 
ordering of the molecules. The orientational elasticity is then described by an energy that depends on Vn 
in a very nonlinear way [CLOO]. There are two key challenges associated with using Vn. 
First, the director is constrained to be a unit vector field because it represents directions, i.e. |n| = 1. 
Second, the head-tail symmetry of nonpolar liquid crystal molecules imposes the constraint that the 
energy is invariant under the transformation n + —n. 


These constraints make the model severely nonlinear. In addition, the constraints — in combination with 
boundary conditions — can lead to the formation of topological defects. Fig. [I|shows classical examples. 
The +1 disclination forms when the boundary conditions set the director on the boundaries. At some 
point in the interior, the director field must be discontinuous if it satisfies the unit constraint; in this 
ideal case, the disclination is at the center by symmetry. There is no way to set up the director field 
to be continuous everywhere with the given boundary conditions. The +3 disclination forms when the 
boundary conditions are as shown in the figure. In addition to the discontinuity at the center, the head-tail 
symmetry also leads to a discontinuous line (in 2D) if it is represented by a director vector field n. If we 
assume that the vectors on the right edge all point upwards and follow the vectors anti-clockwise in the 


upper half and clockwise in the lower half, we find that the vectors point in opposite directions on either 
side of the line (2, < 0,22 = 0). 
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Figure 1: Left: A +1 disclination with a discontinuous director field at the center. Right: A +5 
disclination with a discontinuous director field at the center, and a change in sense along any line 
running out from the disclination. 


The topological nature of disclinations and their connection to boundary conditions is in contrast to in- 
terfaces, cracks, and other defects that are not topological. In the latter, boundary conditions certainly 
play an important role, but through the energy rather than geometry; 1.e., whether a defect forms or not 
depends on the energetics of the particular model. An important consequence of this distinction is that 
phase-field models are readily constructed for the latter class of defects, because adding gradient contri- 
butions to the energy (accounting for surface energies) smooths out the discontinuities. In contrast, there 
is no such simple approach for disclinations which are driven by the boundary. Adding gradient contri- 
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butions to the energy cannot smooth out the defect; on the contrary, the presence of higher derivatives 
makes the energy more singular at a disclination. 


A proposed approach for disclinations is by Ericksen [[Eri91]], in which he introduces a new scalar order 
parameter in addition to the director field. The scalar order parameter is used to track the disclination 
cores and evolves using a separate equation. Another approach in this spirit of introducing another field 
to track defects is [PAD15], where a tensor order parameter is introduced, and a geometrically- 
motivated evolution equation is posed. In contrast to these, we regularize the energy, but do not introduce 
a new field variable nor an associated evolution equation. 


The most widely-applied regularized model that is used to model disclinations is the Q-tensor model 
[IMN14]. Briefly, the order parameter is a tensor Q that is related to the tensor dyad n ® nf} |away from a 
defect, and at a defect the material is no longer uniaxial and cannot be described by a single director. At 
the defect core, the eigenvectors of Q aim to capture some key information about the orientational distri- 
bution function [MN14]. Core effects are introduced by energetic contributions in terms of the gradients 
of Q. Q-tensor models, however, also make certain predictions that are significantly at odds with the 
director-based models. In the static setting, consider 2 infinite parallel plates with tangential anchoring 
and normal anchoring on the top and bottom respectively. As observed by [PMGK94], there are regimes 
in which the Q-tensor model predicts a loss of uniaxiality throughout the domain, i.e., the entire domain 
is defect-like. In contrast, the simple 1-constant director-based model has a solution with a linear vari- 
ation of the angle from top to bottom. Similarly contrasting predictions are observed with the Q-tensor 
model in various geometries and boundary conditions, e.g. [BVD04]. In the dy- 
namic setting, Ericksen and Leslie provided a clear physical model based on angular momentum balance 
for the evolution of the director [Ste04], whereas there is no similar clear principle for the evolution of Q 
in the Q-tensor model. Ideas from kinetic theory are sometimes used, e.g. IYFMWO09], but 
the issue of the closure law for the moments of the orientation distribution function is not clear. 


The issues above motivate the current work. We propose a nonlocal (integral) regularization of the energy 
that accounts for core effects in a natural way, yet retains the many attractive features of a director field 
based model. First, the regularization introduces a length scale that enables us to describe disclinations 
with a finite width; second, using integrals rather than derivatives enables us to have disclinations of 
finite energy even when the director field is topologically constrained to be singular; third, we do not 
need to relax the unit vector constraint as some other models require; and finally, retaining the director 
field description enables the use of angular momentum principles to prescribe the evolution. 


We note that it is particularly important not to relax the unit vector constraint; it is based on the direction 
corresponding to the orientation of the molecules, and hence relaxing the constraint to allow for non- 
unit vectors is not physically meaningful. Therefore, a regularization based on relaxing the unit-vector 
constraint is not particularly physical, even if it is computationally convenient. 


We emphasize 3 important features of the proposed integral regularization of the energy: 


e For homogeneous — i.e. linearly varying — director fields, we should recover the classical energy; 


e Only aconstant director field must have zero energy; any other spatial variation of the director field 
must cost energy; 


e It should not require continuity of the field 


'This ensures the head-tail symmetry is captured. 
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The first requirement ensures that we recover the standard Ericksen-Leslie model when we are far from 
disclinations or in disclination-free configurations. The second requirement ensures that our solutions are 
not “polluted” by unphysical zero-energy soft modes; except for rigid motions, the energy must increase 
when deformed. The third requirement ensures that disclinations can be modeled. 


Our integral model is motivated largely by the peridynamic model of fracture, in which integral operators 
are used to approximate standard elasticity operators, thereby allowing fields to be discontinuous as is 


essential for fracture [|Sil00\|SL10). 


We have ignored the coupling to fluid flow in this work. It is conceptually straightforward to replace the 
classical energy by our regularized energy to construct a model with director evolution coupled to flow. 


1.1 Organization 


The paper is organized as follows. 


e In Section 2] we briefly outline the key features of the energetics of the classical model of liquid 
crystals. This is primarily to show the analogies and contrasts with our proposed model. 


e In Section 3} we describe the proposed nonlocal regularized energy without gradients. Specifically, 
we (1) motivate the connection between integral operators and gradient operators, (2) formulate the 
nonlocal analog of the classical energy, (3) discuss the evolution equation and boundary conditions, 
and (4) extend the energy to account for head-to-tail symmetry. 


e In Section |4| we incorporate the effects of different moduli for bend and splay. 


e In Section |5} we examine the structure of disclinations, and the role of bend and splay in the 
stability of integer vs. half-integer disclinations. 


e In Section|6| we examine the coalescence of a +3 disclination and —} disclination coming together. 


1.2. Notation, Definitions and Values of Model Parameters 


Boldface denotes vectors and tensors. We have used Einstein convention, i.e. repeated indices imply 
summation over those indices, except when noted. 


Hats appended to functions refer to their value at a, g, ..., Le. for a field f(a), we write f = f(x), f = 
ice ee 


1 = a3 — Q2 => 0 is the rotational or twist viscosity of the director. It sets the relaxation time of the 
director field. We set 7; = 1 in our model system. a3 and az are the Leslie viscosities [|Ste04]. 


6 is a characteristic defect lengthscale, that we set to 5. Numerical calculations use a grid spacing of at 
least 1. 


Cs(\a — &|) is the nonlocal weight function; we assume the form 


1 . “a 
_an— J -@ ifle-—a| <6 


4 
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The normalization is appropriate for 2D. 


2 Classical Model of Nematic Liquid Crystals 


We briefly summarize the key features of classical nematic director dynamics to enable us to point out 
the key differences and similarities with our proposed nonlocal model. The classical model of nematic 
liquid crystals was formulated by Frank, and Ericksen and Leslie, and is very well described in [[Ste04]. 
For simplicity, we neglect in this paper the effect of flow, i.e., we assume that the only dynamics is due 
to director reorientation. In addition, we use the standard assumption that the inertia associated with 
director reorientation is negligible. 


The Frank free a density ( is in general a function of the local values of n and Vn. The energy 
of the body 2 is E [n = fou n, Vn) dQ. We fix n or set €,17;,;217; = 0 on the boundary 0.2, and 
[re] =, 

In the simplified 1-constant approximation, w = £|Vn|?. 

The classical approach relies on angular momentum balance [Ste04] to derive the evolution equation for 
mn. Equivalently, we can derive this using a gradient descent approach. We prefer the latter method here 


because it is easier to adapt to the nonlocal energy because we do not have to first identify a nonlocal 
torque, but can instead do so from the final form of the equation. 


The gradient dynamics has the form 7,(7,6n) = — £E [n+ ebn]|__ >: To preserve |n| = 1, the 
variations in m must have the form 6n; = €:;£7;D,%, Where p; is arbitrary. The dynamical equation is then: 
Ow Ow Ow Ow 
wxn=— | — —div——] xno y7n=- — di + 2.1 
ae (+ _ a oe (5 ad aa) i 


where \ is the Lagrange multiplier conjugate to the constraint |n| = 1. Enforcing n-n =1>n-n=0 
to eliminate \ gives the dynamical equation: 


dw. dw \ | Ow _ Ow 
wn =— (5 — div a (n. ah —n-div wa) n (2.2) 


The dissipation associated with this dynamics is: 


dE Ow Ow 
a di -n dQ 
dt fe ($= ~ a) 
1 Ow dw \/? Ow dw \? 
_ _ dj = — —_n-di > 
(= div a) (n. n- div ) dQ >0 


an on ovn 
The dissipation is strictly non-negative?| 


(2.3) 


V1 Jxen 


Equilibrium configurations are obtained by setting m = 0. In the 1-constant energy, this gives €;,174,;;71 = 
0. 


> Gradient dynamics with a unit vector constraint is not necessarily dissipative. The Landau-Lifshitz-Gilbert (LLG) dy- 
namics in micromagnetism (e.g., [AKST14]) is given by (n,p) = — 4E [In + cdn] = with dn; = jr; Pk, Where pz, is 
arbitrary, in the undamped case. The inner product above is not the same as the Ericksen-Leslie dynamics, and it can be shown 
that the energy in the undamped LLG dynamics is conserved. 
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Alternatively, we can parametrize n = (sin 6(a), cos @(a)) which automatically satisfies the unit vector 
constrainf}} Then, ni, = cos(@)0; and no, = —sin(@)6,;, giving for 1-constant energy: 


E{9] = > [ 0.0; d2 (2.4) 


The equilibrium configuration is described by @;,; = 0 in (2, with boundary conditions 6 = 6 or 6,7; = 0 
at each point on the boundary 02. The description in terms of @ has the advantage that the unit vector 
constraint is exactly satisfied at all times for numerical solutions that discretize in time and use finite time 
steps. The description in terms of n will require a careful projection scheme at each time step to bring 
the director back to the unit sphere. Our approach for the nonlocal energy will be to formulate in terms 
of a unit vector field n, but use @ for the numerical solution. 


The head-tail symmetry is important to describe half-integer disclinations: even away from the core, 
there is a line (in 2D) across which the sense changes by 7. Using either @ or a vector field for n does 
not allow us to do this easily with the classical energy, but will be possible with the proposed nonlocal 
regularization. 


3 Nonlocal Regularized Free Energy Without Gradients 


Our nonlocal regularized energy is motivated by the peridynamic model for fracture with 
some key differences. Peridynamics has provided an important model for fracture calculations because 
of the key feature that spatial derivatives are not required in the theory, hence allowing for discontinuous 
displacement fields. 


3.1 A Motivating One-Dimensional Example 


Consider the operator D that acts on a function f as defined: 


x+6 
Bfla) = 55 | (Fe) - F@)(@-) ds G.I 
When f(x) = Ax + B, we evaluate Df(x) = A. For linear functions, Df matches the derivative. 
Consider the Taylor expansion f(x) — f(#) = (# — 4) f'(x) + $(a — #)?f" (x) + .... Then, Df (2) = 
fi(x) + C362,f’"(x) + ..., where the even derivatives vanish by symmetry. Therefore, the operator D 
provides an approximation to the derivative when 6 and the higher derivatives are sufficiently small. 
However, when /f is not smooth such as at a discontinuity or other singularity where the derivatives are 
not well-defined, the original integral expression continues to be well-defined. 


Physically, the operator can be understood as averaging the slope of f over a neighborhood of x of size 
0, 1.e., averaging fe Fe) However, to avoid singularities from the denominator, we use the weight 
(2 — #)? in the process of averaging. When the slope is constant, i.e. a linear function, we recover exactly 
the derivative; when the function is slowly varying, we get higher-order regularizing terms; and when the 
function does not have well-defined derivatives, we obtain a finite value that is physically related to the 


average slope in the neighborhood. 


3This is in 2D, but the extension to 3D is straightforward. 
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While Df provides an attractive alternative to the derivative df /dx, there is one key drawback. The 
derivative df/dx = 0 only when f = const. and not otherwise. However, while Df = 0 when f = 
const., it can also be 0 for other functions; in fact, any non-constant function that is antisymmetric about 
x will have the slopes exactly cancel out to give an average value of 0. Therefore, one can have spurious 
soft modes with zero energy that have no physical basis. 


To get around this, we first observe that the free energy density consists of terms that are all components 
of the gradient raised to the second power [|Ste04]. In the specific case of the 1-constant energy, the 
energy density is simply the gradient squared. In one dimension, this would require us to approximate 


2 : . —_— : 
(£) . Motivated by the discussion immediately above, the operator 


Dof (x) = / . (fc f8y’ di (3.2) 


— t—-DL 


would provide a nonlocal regularized approximation up to normalization. The integrand is always non- 
negative, thereby not allowing any cancellations. For D2 f(a) to evaluate to 0, the integrand must be 0 
everywhere, ensuring that / = const. is the only possibility. 


3.2 Nonlocal Regularization in Three Dimensions 


Following ideas from peridynamics [|Sil00], we write the nonlocal analog to the Frank energy as: 


_K Fs cas 
=f | OG aw a ae (3.3) 
xen Jaen 2 |jz — 2| |e — x| 


where C’; is a function that is symmetric in the arguments. 


What motivates this expression? The classical gradient is the vector pointing in the direction in which 
the increase of 7) is maximum. Our expression examines the change in n between & and every point z in 
2, and weights the direction « — x by the change in that direction; the larger the increase, the greater the 
weight for that direction. So the direction of largest increase features most prominently in the average. 


It is natural for points # that are closer to x to play a larger role in the averaging, and this weighting 
is introduced by the function C;(x, a). This introduces a lengthscale, denoted 6, into our model that is 
not present in the classical Frank free energy. In this paper, we will restrict ourselves to isotropic and 
spatially homogeneous systems, and hence we will write Cs(|a — #|). Further, we use Cs to introduce a 
cut-off distance, i.e. C5(|a — &|) = 0 when |w — a| > 0. Finally, for conciseness, we absorb = =k Z = 
in (3.3) within Cs by redefining it appropriately, but retaining the same symbol for sonvenionce. Putting 
this together, we can write: 


nl=> | 2G (ja —#|)-(m%—n)*? dV~ dV, (3.4) 
LreEQI& 


en 2 


We now consider a physical interpretation of (3.4). From a physical perspective, the classical Frank 
model states that there is an energetic cost to molecules not being aligned with each other. Loosely, 
it penalizes changes in the molecular orientations at “adjacent” material points, and in the limit this 
becomes the gradient. Our model compares the molecular orientation at a point x with the molecular 
orientation at every point z contained in {2. To leading order, the energy contained in the interaction 
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between molecules at these points is $C5(|a — #&|) (n(x) — n(a))”, where Cs(a, &) is analogous to a 
spring constant and depends on the distance between the points. The total free energy simply sums over 
all pairs of molecules in {2, with an additional factor of 5 to correct for double-counting of bonds. While 
molecular interactions are typically very local, in homogeneous non-defect regions we can think of the 
molecular interactions as being renormalized. 


3.3 Accounting for Head-Tail Symmetry 


Typical nematic liquid crystal molecules are non-polar and have head-tail symmetry. That is, chang- 
ing nm «++ —n should not change any physical quantities. This is particularly important in half-integer 
disclinations, e.g. Fig. |1} where there is a plane — or line in 2D — across which the sense of the director 
changes. If represented by a unit vector field, the vector flips head-to-tail; if represented by @ in 2D, there 
is discontinuity of 7. These discontinuities should not contribute to the energy or equilibrium equation. 
However, the energy in is sensitive to such a transformation. We now discuss how to modify the 
energy to reflect the physics. 


Consider two directors that are within 6. The minimum energy configuration is n; - m2 = 1 (they are 
parallel). The highest energy configuration is n; -m 2 = O (they are normal). We further want that the 


configuration 2; - m2: = —1 (they are anti-parallel) is a minimum energy state, with the same energy 
as the parallel state. Some choices of energy that satisfies this are —|n1 - n2| = —|cos(O; — 62)], 
—(n1 +12)? = —cos?(0; — 62), ete. 


Consider now the integrand from (3.4). We can write (n — 7)? = 2—2n-% = 2—2cos(O—6). Since it is 
an energy, we can ignore constants and write the energy as —2n- 7% = —2.cos(9—6). From the discussion 
in the previous paragraph, it is clear that this does not behave well when n and 7 are anti-parallel, and 
shows the problem with @3.4). 


Based on these physical considerations, we replace the integrand in nonlocal energy density by — |2 — (n — n) 


—2|n-n|. That is, we propose the following nonlocal regularized free energy: 


_K 
=i | 1 oy( (| — &|)--|2—(n—7)"| dVe dVo (3.5) 
LEQSE 

SS 


en 2 


free energy density at x 


This recovers the classical model when the anti-parallel nature is not an issue because it is quadratic in 
m—. It is obvious that our proposed expression does not impose any continuity requirements on n. The 
energy density — and therefore the total energy — is bounded and well-defined even if n has all kinds of 
discontinuities, as long as n itself remains bounded; this is automatic due to the unit vector constraint. It 
is also obvious that the energy density is non-zero except for constant fields if we restrict C’s(|a—a|) > 0 
because the integrand is non-negative. 


We could, in principle, use many other expressions that resolve the head-tail issue, e.g. —2(n - n)? for 
the free energy density. But these would not recover the classical model in the limit. 


3.4 Equilibrium and Evolution Equations 


The equilibrium equation is obtained from setting the functional derivative of the energy from (3.5) to 


0. That is, 4F [n+ cbn)|__, = 0. To satisfy the unit vector constraint on n, we use the variation 


Disclinations without Gradients: A Nonlocal Model for Topological Defects in Liquid Crystals R. Macedo, H. Pourmatin, K. Dayal 


ON; = €ijnN;Pk Where p is an arbitrary vector field. This gives: 


K ; . ; _ 
as a i Cs(\a — &|) - sign(2 — (n — ny”) - Eijk (M4 — fi) (NP — MyPr) AV Ve 
rEQ J BED 


K 
=e= / | C;(|a — &|) - sign(2 — (n — f)”) - e:jn (ni — A) NPE UVa UV 
2 2EQ J BEN 


K (3.6) 
dicate | i Cs(\a = £|) ‘ sign(2 = (n cae n)’) * €ijk (n; = Nj) N;Pk dV, dV, 
wEN JEN 
= K | mf Cs(|a — &|) - sign(2 — (n — f)”) - e:jn (ni — 74) Nj Va dV 
xen LEQ 
where we have relabeled x + x to combine integrals in the second line. 
Using that p; is an arbitrary field, we eliminate the integration over x: 
Eijk if Cs(\a — &|) - sign(2 — (n — 2)?) - (nk — Ax) Va = 0 (3.7) 
&EQ 
We compare this with the statement of equilibrium in the classical 1-constant model: €;;,2;N%,11 = 0. 
In 2D, an orientational description n = - 0(x), cos O(a)) gives: 
_*k ‘ 
=i | Lot (ja — 2]) (1 — | cos(? — 6)!) dVz, dV (3.8) 
we JLEN 2 
for the energy, and equilibrium configurations are described by: 
K | Cs(|x — &|) sign [cos(4 — 6)| sin(9 — 0) dVz = 0 3.9) 
&EQ 
As in the classical case, we deduce the evolution equation using y;(n,dn) = — LE [n+ €dn] = using 
the variation dn; = €j.;Px, Where p is arbitrary. This gives us: 
ynxn= (« | C;s(|a — &|) -sign(2— (n —f)*)-(n—") dV WW) xn 
BEN (3.10) 
=n. = K | Cs(|a — &|) - sign(2 — (n — n)”) -(n — 1) dVg dVz + An 
BEQ 
Enforcing n-n = 1 > n-n = 0 to eliminate X gives: 
yn = U-nen)-K f Cs(|a — &|) - sign(2 — (n — ”)*) - (n — 1) dVe dVx (3.11) 
BEL 
The corresponding 2D orientational description is: 
_ K | C5(|a — #|) sign [cos(6 — 6)| sin(9 — 0) dVz = 0 (3.12) 
BEL 
The dissipation associated with the nonlocal model can be computed: 
-@ = Kk ; 
=i | C;s(|x — &|) - sign (2 — (n — ”)”) -(n—”)- (n—7) dV; dV, 
wE2Q JZEN 
ef fe O;(|a — &|) - sign (2 — (n — n)”) + (n—n) dVz dV, (3.13) 


=e ae (I-n@®n)-Adv, 
Va 
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where A := K J.,._,C3(|a — @|) - sign (2—(n— n)”) -(n—7) dVz. Since I — n & n is positive 
semi-definite, it follows that the dissipation is always non-negative. 


3.5 Boundary Conditions 


As in the peridynamic model of elasticity, boundary conditions cannot be applied in the classical sense. 
Roughly, imposing conditions on the boundary — a set of measure 0) in the integral — provides an in- 
finitesimal contribution to the integral. However, a physically-natural solution to this issue is to instead 
impose conditions over a layer of finite thickness — with thickness of order 6 — on the boundary. Given 
that in our model molecules interact over a finite range, it is natural to provide boundary conditions with 
more structure than holding them fixed on a low-dimensional set (the boundary). For instance, to impose 
that the nematic molecules have a specific orientation on the boundary, we hold the orientation fixed for 
all molecules within the boundary layer region. That is, O(a) = 09 Va € (Qy where {2,; is a portion 
of the §2 with finite volume. In this paper, this is the only type of boundary condition that we use, but 
other boundary conditions can be similarly smeared-out following the techniques used in peridynamics 


S110} |SL10) |DBO6]. 


4 Bend and Splay: Extension Beyond the 1-Constant Energy 


In the classical 1-constant model, a single modulus Kk defines the energy in the bend, splay, and twist 
modes. The more general Frank free energy has different moduli for each of these modes. We examine the 
same issue in the nonlocal regularized energy. We restrict ourselves to 2D and therefore only distinguish 
between bend and splay, but an extension to 3D that also incorporates twist is conceptually analogous. 


The physical basis for 3 different moduli is seen from Fig.|2} We can differentiate between these modes as 
follows. Splay can be differentiated from both bend and twist by noticing that the quantity (n—7)-(a—z) 
is 0 for both bend and twist, but is nonzero — and can be normalized to 1 — for splay. Similarly, bend can 
be differentiated from both splay and twist by using the quantity n - (a — x), which is 0 for both splay 
and twist but is nonzero and normalizable to 1 for bend. 


Since we are currently only working in 2D and considering only bend and splay modes, we need only 
differentiate these 2 modes here. We write for the modulus kK = Kk, (1 — In : Six) + Ky In : xl, 


= =— aE the coefficient for splay and Kx; is the coefficient for bending. When n - dx = 0, 


then we have ane and when n - 6a = 1 we have bending. 


Using the orientational description for compactness, the energy can now be written: 


EI) = [ RO.@:8) | o(|e — al) (1—| cos(@ — 8))) dVz dV (4.1) 
wEQ JEN 
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Figure 2: Bend, Splay and Twist modes. These modes can be differentiated 
in our nonlocal model by using the quantities (n — n) - (a — &) and n 


(a — x). When normalized, these quantities are both 0 for twist, 0 and 1 
respectively for bend, and 1 and 0 respectively for splay 


We can also sid the bending and splay contributions individually 


swim ff 
LEQSE 


(|x — 21) In . bax (1 — | cos(6 — 0)!) dVs dV 
en 2 


bend (4.2) 
+ sk, Lbs 5 O5( |a — x) )(a- In - dar) (1- |cos(6 — 6)]) dVz, dV x 


splay 


The evolution equation is obtained from a gradient descent based on 6 without any constraints 


_ [calle — £|) (FS (1 —|cos(6 — 0)|) + 


K(@) - sign [cos( 4 — 6) 
The equilibrium equation is obtained by setting the right side above to 0 


sin( — ®)) dVz = 0 
(4.3) 
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5 Structure and Stability of Individual Disclinations 


We use the model to examine the stability and structure of various half-integer and integer disclinations. 
We examine these using the 1-constant nonlocal energy as well as the energy that differentiates between 
bend and splay. 


5.1. The 1-Constant Nonlocal Energy 


In the 1-constant nonlocal energy, we find that only +4 disclinations are stable; disclinations of higher 
strength always split into the appropriate number of +3 disclinations. We show some examples of this 
below. 


In all of these calculations, we use a discretization mesh size of 1 with 6 = 5, except where stated 
otherwise. We use a square domain with size 200 x 200. For the starting configuration, we use the 
classical formula for a disclination n(x) = (sin(k0o(a) + a), cos(ko(a) + a), where tan 09 = r2/21, 
and then evolve using the equations presented above. The initial and final configurations are shown in 


Fig. 


This finding agrees qualitatively with calculations in [MHML05] |BKZ98] which use the Q-tensor and 
molecular dynamics methods respectively, in that only +4 defects are stable in the 1-constant energy. It 
is consistent also with the heuristic that the energy of a disclination scales with the square of the charge 


that it carries [dP95], as well as rigorous work [|BPP12]. 


While the final states are largely in agreement with other aproaches, the dynamics is qualitatively different 
in our model compared to the Q-tensor calculations from [MHMLO5]. They find that a +2 defect first 
splits into two +1 defects, and then splits a second time into four 15 defects. In our calculations, all of 
the defects split almost instantly into +4 disclinations that are then repelled from each other and move 
apart. While their geometry is circular and hence different from our square domain, in both their work 
and ours the domain is large enough that it is unlikely that this is the reason for the difference. It is likely 
a consequence of the differences between the evolution equations in the Ericksen-Leslie model with our 
energy vs. the Q-tensor model. 


We also notice that the energy (3.4) has stable +1 disclinations and does not form +3 disclinations. 
However, this occurs due to the unphysical feature that the lack of head-tail symmetry associates a very 
high energy to the line (in 2D) across which the director changes sense that accompanies +3 disclinations. 


We then examine the detailed structure of the core for +4 defects. For these calculations, we use a 
much finer discretization of 0.1 with the same value of 6. In an initial set of calculations, we began 
with the classical solution as above, and evolved to find the relaxed configuration. The solutions in our 
model are extremely close to the starting classical solution, and the changes are dominated by numerical 
discretization errors. So we instead use as initial the classical solution everywhere outside a core region of 
radius 2; within this region, we set all directors to point vertically. The director field relaxes to a solution 
that is again very similar to the classical solution, but the disclination translates during the relaxation due 
to the severe perturbation that we have induced in the initial conditions. We therefore perform a best fit of 
the relaxed configuration against the formula for a translated disclination, treating the disclination center 
as the parameters to optimize over. Fig. |4]shows the change in the director angle. The difference between 
our solution and the classical disclination solution is confined to the core region, and even there is small. 
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Energy 


Figure 3: Top: A +2 disclination constructed using the classical solution (right), and the final 
configuration of four well-separated +5 disclinations (right). Middle: The similar process for a 
+3/2 disclination decomposing into three +3 disclinations. Bottom: A —1 disclination decom- 
posing into two —s disclinations. The vector field is overlaid on the energy density. We observe 
similar decompositions in +1, —3/2, and other disclinations stronger than +3. 
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30 40 50 60 70 80 
x1 


Figure 4: Angular difference between a classical disclination solution and our relaxed disclination 
configuration for +5 (left) and —4 (right). 


Change in theta (radians) 


5.2 Disclinations in the Bend-Splay Energy 


In the previous subsection, we see that the nonlocal 1-constant energy has only +3 disclinations. In 
experiments, however, it has been noticed that +1 defects dissociate into two closely-spaced +3 cores 
[MPR14]. As described in [YRO2], the dissociation of +1 disclinations is controlled by the ratio of bend 
to splay energies. Looking at Fig. || we see that the +1 disclination is dominated by bend, whereas 
the +5 disclination has a significant amount of both splay and bend. We notice that in our approach to 
differentiating between bend and splay based on n - (a — a), there will be some level of splay character 
in the classical +1 defect because of contributions from points that are at around +7//4 to the director 
orientation. We examine the behavior of a +1 disclination in our model at two limits, one with kK, = 
10K, and the other with K, = 0.1K>. 


We first plot the partitioning of energy between bend and splay in the classical solution with kK, = K, in 
Fig. |5| As we expect, the bending energy is generally dominant, though the splay energy is not vanishing 
as explained above. 


We then evolve the director field from this initial configuration. In the case when the bending modulus 
is large compared to the splay modulus (KX, >> K,), bending is expensive energetically and hence the 
system dissociated into two 15 disclinations in which some of the distortion can be accommodated as 
splay (Fig. [6] top). In the case when the bending modulus is small compared to the splay modulus (kK, < 
K,), bending is not expensive and hence the system retains roughly the original disclination configuration 
thereby reducing the amount of splay (Fig. [6] bottom). Intermediate configurations between these two 
limits lead to partially dissociated disclinations that are qualitatively comparable to the observations in 
and the theoretical model in [YRO2]. 


14 


R. Macedo, H. Pourmatin, K. Dayal 


Disclinations without Gradients: A Nonlocal Model for Topological Defects in Liquid Crystals 


ABsau3 


bid 


ith kK, = K,. The 
lination core. 


1eS W1 


to splay (left) and bend (right) energ 


ficantly higher energy. Both are local 


in 


f the energy 


Partitioning o 
distortion has s 


Figure 5 


1SC 


d around the d 


1Ze 


igni 


bending 


A619U3 


EEX X XNA 

PAPA A A EEE 

LAAA A Anne, ERR RANA AY 

GLEE ARAN AANA 
er oe 

Af AA Ao. Ss 


NAAAAAS Se 
SANA 


SSS 
\NSSS SS 


into splay (left) and bend (right) in the final configuration 


Energy partitioning 


: Top row: 


Figure 6 


for K, = 1044, i.e., splay is much more expensive than bending. The defect retains a +1 character 


in 


lay (left) and bend (right) 


ion. Bottom row: Energy partitioning into sp 


1stort 


lay d 


to minimize sp 


han splay. The defect 


ive t 
her splay d 


h more expens 


ing is muc 
d compensates with ah 


the final configuration for K, = 10K4, i.e, bend 


istortion. 


1g 


istortion an 


ing d 


issociates to reduce the bendi 


d 


15 


Disclinations without Gradients: A Nonlocal Model for Topological Defects in Liquid Crystals R. Macedo, H. Pourmatin, K. Dayal 


6 Coalescence of 1s and 5 Disclinations 


We next use our model to examine the evolution of 2 disclinations of equal-and-opposite +4 topological 
charge. Two such disclinations will attract each other, causing them to move towards each other and 
eventually coalesce to form a disclination-free material. Whether these defects move at the same velocity 
towards each other, and hence coalesce at the midpoint between the initial defect positions, has been an 
important question. This has been the focus of both experimental and theoretical investigation. A key 
experiment observing coalescence of point defects in nematics shows that there is asymmetry in their 
motion [CB03]. Closely related situations have been studied numerically 'SZ02] and 
theoretically [GSV02]; while some of these systems are significantly different from the nematics studied 
here, a key broad finding is the significant influence from hydrodynamic effects. 


A priori, there is no reason for equal-and-opposite topological disclinations to move at the same velocity. 
Disclinations are not analogous to electrical point charges which have no internal structure and are there- 
fore perfectly anti-symmetric. Rather, the topological charge is simply one average measure — admittedly 
an extremely important measure — of the defect structure. The detailed structure of the disclination plays 
an important role in the dynamics, and even if different disclinations of the same strength have the same 
force acting on them, the detailed structure is essential in setting the dynamics in response to the force. 
Returning to the electrical analogy, a complex charge distribution can be characterized by the net charge, 
but the dynamics of the charge distribution in response to another charge distribution cannot be charac- 
terized solely in terms of the net charges of both distributions. 


We examine this situation using our simple |-constant regularized model using ++ disclinations. Fig. 
[7| shows the initial configuration with two +4 disclinations and the final defect-free configuration. The 
force between disclinations falls off very rapidly with distance, and hence there is initially little move- 
ment. However, as the disclinations begin to come closer, the force increases, thereby further increasing 
the velocity, and so on. Eventually, the disclinations coalesce and leave behind a disclination-free ma- 
terial. A movie of the coalescence process is part of the supplementary material, and we quantify the 
approach in Fig. While small, there is an unmistakable difference in the velocities of the positive 
and negative defects. Given that our l-constant regularized model contains no complexities such as 
flow, or even different moduli for bend and splay, it suggests that this asymmetry in disclination veloc- 
ity has its origin in the asymmetry of the director field; flow further enhances this effect when present 


(TDY02)|SZ02). 


7 Conclusions 


We present a nonlocal regularized model that exploits integral operators to achieve a generalization of the 
classical Ericksen-Leslie model for nematics. Near disclinations, the model regularizes defect cores, and 
as the distortion becomes uniform, it recovers the Ericksen-Leslie model. The use of integral operators 
—as opposed to differential operators in the classical approach — enables us to model situations in which 
the disclination field is not continuous. As we describe, the integral operators tend to the differential 
operators in a physically meaningful way. 


Our approach uses in an essential way the important ideas behind peridynamics [|Sil00} |SL10]. How- 
ever, the head-tail asymmetry that is essential in describing liquid crystals does not have an analog in 
peridynamics which is tailored to deformation fields. For this reason, it is not possible to directly use 
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Figure 7: Initial configuration (left) with two + 
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generalizations of peridynamics to construct energies in liquid crystals. In particular, the 3-constant 
Frank free energy that we construct appears quite different from, and does not seem to have analogies to, 
strategies used in peridynamics to model complex materials. 


This is a preliminary step towards the formulation of this class of models. While we have demonstrated it 
in the context of nematics, it is readily generalizable to other liquid crystalline material such as smectics. 
The key issues that we have dealt with here — namely unit vector fields that have discontinuities — are the 
central challenges in this class of materials. In addition, we have not fit the many material parameters 
here to specific materials. For instance, the regularization scale 0, the precise form of the kernel C's, and 
the precise angular dependence of /, can all be tailored to capture specific types of physics; for instance, 
if detailed dynamic behavior is to be predicted. We have also restricted ourselves to 2D for simplicity, 
but an extension to 3D is conceptually simple. Such an extension would enable the interrogation of more 
realistic geometries, as well as enable us to examine twist distortions. 


Our numerical approximations have used finite differences for simplicity. However, powerful numerical 
techniques such as discontinuous Galerkin finite elements (DGFEM) can enable the efficiency of adaptive 
meshing and finite elements to attack much bigger problems and complex geometries with our model. 
While defect tracking can be enormously expensive, results from suggest that DGFEM 
offers an efficient approach: it is a conforming approximation for peridynamics, and the results appear 
not very sensitive to the errors in tracking the defect. 
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